BSLE_P_Code
Morgan Barnes
‘2023-08-18’r Sys.Date()’
##This is the code for figures and statistics from Barnes et al. (2024): phosphorus composition along a burn severity gradient ##Contact: Morgan E Barnes morgan.barnes@pnnl.gov ##Make your environemnt
Set Up
Load libraries
require(pacman)## Loading required package: pacmanpacman::p_load(tidyverse,
dplyr,plyr,ggplot2,
lubridate, tidyr, tidyverse,
scales, cowplot, gridExtra, RColorBrewer,
wesanderson,readxl,openxlsx, stringr, hablar,
janitor, stringr, naniar, forstringr, reshape2,
emmeans, rstatix, ggpubr, Hmisc, corrplot, ggrepel,
stringr, lme4, lmerTest, car, performance, vegan, psych, lavaan)Confirm working directory:
## [1] "/Users/barn772/OneDrive - PNNL/Documents/GitHub/rcsfa-RC3-BSLE_P/scripts"##Load data The datasets and metadata needed for P BSLE Gradient manuscript are: 1) Chemistry datasets solids and liquids w/ all metadata that can be used directly from data package: This includes: a) Total elemental composition of solids b) Total elemental composition of leachate particulate and filtered fractions c) Molybdate reactive P of leachate filtered fraction d) XANES of solids e) NMR of solids
- Chemistry data that needs post-data package processing after download from BSLE Data Package v3; this was done and provided within the BSLE P Manuscript Data Package:
- XANES of solids (spectra of samples and reference compounds and linear combination fits)
- NMR of solids (deconvolution)
Download the BSLE Data Package
Until we can get the API to work, do this manually and store on your
local. I recommend storing inside your github folder and gitignoring the
subfolder with the data in it by adding this to your .gitignore:
#data folder data/
BSLE v3 data package URL: https://data.ess-dive.lbl.gov/view/doi:10.15485/1894135 BSLE P manuscript data package URL: INSERT HERE!!! fix
Bring in the data package data:
#metadata
Meta <- read_csv("../data/BSLE_Data_Package_v3/v1_BSLE_Metadata_and_Protocols/BSLE_Burn_and_Laboratory_Metadata.csv") %>%
dplyr::rename(Sample_ID = Sample_Name,
Leachate_vol_mL = Leachate_Volume,
Char_leached_g = Char_Leached) %>%
naniar::replace_with_na_all(condition = ~.x ==-9999) %>%
naniar::replace_with_na_all(condition = ~.x == "N/A") %>%
mutate(Burn_Severity=replace(Burn_Severity, Burn_Severity=='unburned', 'Unburned'))## Rows: 234 Columns: 26
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (19): Burn_Code, Parent_ID, Sample_Name, Feedstock_Species, Land_Coverag...
## dbl (7): Max_Burn_Temp, Dry_Fuel_Before_Burn_Weight, Approx_Char_After_Burn...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Meta<-Meta[ c(2,3,5,8,21,25,26,6,20,22,23)]
Meta$Burn_Severity=factor(Meta$Burn_Severity, levels= c("Unburned", "Low", "Moderate", "High"))
#chemical data from BSLE v3 data package
ICP_Solid <- read_csv("../data/BSLE_Data_Package_v3/v3_BSLE_Data/BSLE_ICP_Solid.csv")## New names:
## Rows: 33 Columns: 12
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (12): #Columns, 11, ...3, ...4, ...5, ...6, ...7, ...8, ...9, ...10, ......
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...3`
## • `` -> `...4`
## • `` -> `...5`
## • `` -> `...6`
## • `` -> `...7`
## • `` -> `...8`
## • `` -> `...9`
## • `` -> `...10`
## • `` -> `...11`
## • `` -> `...12`ICP_Leachate <- read_csv("../data/BSLE_Data_Package_v3/v3_BSLE_Data/BSLE_ICP.csv")## New names:
## Rows: 128 Columns: 12
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (12): #Columns, 11, ...3, ...4, ...5, ...6, ...7, ...8, ...9, ...10, ......
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...3`
## • `` -> `...4`
## • `` -> `...5`
## • `` -> `...6`
## • `` -> `...7`
## • `` -> `...8`
## • `` -> `...9`
## • `` -> `...10`
## • `` -> `...11`
## • `` -> `...12`MBD <- read_csv("../data/BSLE_Data_Package_v3/v3_BSLE_Data/BSLE_MBD.csv")## New names:
## Rows: 71 Columns: 5
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (5): #Columns, 4, ...3, ...4, ...5
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...3`
## • `` -> `...4`
## • `` -> `...5`ICP_PNMR <- read_csv("../data/BSLE_Data_Package_v3/v3_BSLE_Data/BSLE_ICP_PNMR.csv")## New names:
## Rows: 33 Columns: 20
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (20): #Columns, 19, ...3, ...4, ...5, ...6, ...7, ...8, ...9, ...10, ......
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...3`
## • `` -> `...4`
## • `` -> `...5`
## • `` -> `...6`
## • `` -> `...7`
## • `` -> `...8`
## • `` -> `...9`
## • `` -> `...10`
## • `` -> `...11`
## • `` -> `...12`
## • `` -> `...13`
## • `` -> `...14`
## • `` -> `...15`
## • `` -> `...16`
## • `` -> `...17`
## • `` -> `...18`
## • `` -> `...19`
## • `` -> `...20`npoctdn <- read_csv("../data/BSLE_Data_Package_v3/v3_BSLE_Data/v2_BSLE_NPOC_TN.csv")## New names:
## Rows: 328 Columns: 6
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (6): #Columns, 5, ...3, ...4, ...5, ...6
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...3`
## • `` -> `...4`
## • `` -> `...5`
## • `` -> `...6`SolidCN <- read_csv("../data/BSLE_Data_Package_v3/v3_BSLE_Data/BSLE_CN.csv")## New names:
## Rows: 83 Columns: 6
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (6): #Columns, 5, ...3, ...4, ...5, ...6
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...3`
## • `` -> `...4`
## • `` -> `...5`
## • `` -> `...6`pH <- read_csv("../data/BSLE_Data_Package_v3/v3_BSLE_Data/BSLE_pH.csv")## New names:
## Rows: 247 Columns: 5
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (5): #Columns, 4, ...3, ...4, ...5
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...3`
## • `` -> `...4`
## • `` -> `...5`NMR_spectra <- read_csv("../data/BSLE_Data_Package_v3/v3_BSLE_Data/BSLE_31P-NMR.csv")## Rows: 374667 Columns: 58
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (58): Chemical_Shift_parts_per_million, BSLE_0001-solid, BSLE_0002-solid...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.#chemical data from BSLE P manuscript data package
NMR <- read.csv("../data/Barnes_2023_BSLE_P_Gradient_Manuscript_Data_Package/NMR_Deconvolutions.csv") %>%
filter(Degradation_Corrected == 'Yes')
XANES <- read.csv("../data/Barnes_2023_BSLE_P_Gradient_Manuscript_Data_Package/P_XANES_Samples/P-XANES_LCF.csv")
XANES_spectra <- read_csv("../data/Barnes_2023_BSLE_P_Gradient_Manuscript_Data_Package/P_XANES_Samples/xanes_example_spectra.csv") #this file was generated using XANES_Spectra_Merging script ## New names:
## Rows: 326 Columns: 9
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," dbl
## (9): ...1, eV, BSLE_0013-solid_P-XANES.xmu.nor, BSLE_0007-solid_P-XANES....
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...1`Clean up data functions:
essdive_template_cleaning <- function(dataframe){
cleaned <- dataframe %>%
janitor::row_to_names(row_number = which(dataframe[,1] == "Field_Name")) %>% #this assumes that this name is somewhere in column #1. That follows some (most?) ESS-DIVE templates.
slice(which(Field_Name == "#Start_Data"):which(Field_Name == "#End_Data")) %>%
select(!Field_Name)
}
noreporting_format_cleaning <- function(dataframe){
cleaned <- dataframe %>%
janitor::row_to_names(row_number = which(dataframe[,1] == "Date")) #this assumes that this name is somewhere in column #1. This works for the thermocouple datasets
}pH_clean <- essdive_template_cleaning(pH) %>%
dplyr::rename(Sample_ID = Sample_Name,
pH = '00403_pH') %>%
mutate_at(c(3), as.numeric) %>%
separate(col=Sample_ID, into = c("Sample_ID", "FilterSize"), sep="-") %>%
filter(grepl('unfilt', FilterSize)) %>%
select(Sample_ID, pH)
#wrote modifies so -9999 is NA, below LOD is NA, and if above or below standard curve then I pull out those values
npoctdn_clean <- essdive_template_cleaning(npoctdn) %>%
dplyr::rename(NPOC_mg_C_per_L = `00681_NPOC_mg_per_L_as_C`,
TN_mg_N_per_L = `00602_TN_mg_per_L_as_N`,
Sample_ID = Sample_Name) %>%
select(Sample_ID, NPOC_mg_C_per_L, TN_mg_N_per_L) %>%
filter(grepl('filt0.7', Sample_ID)) %>%
mutate(Sample_ID = stringr::str_remove(Sample_ID ,"-filt0.2"),
NPOC_mg_C_per_L = case_when(NPOC_mg_C_per_L == -9999 ~ NA,
str_detect(NPOC_mg_C_per_L, "LOD") ~ NA,
str_detect(NPOC_mg_C_per_L, "Below_0.1_ppm") ~
str_extract_part(NPOC_mg_C_per_L, "_ppm_Below_0.1_ppm", before = TRUE),
TRUE ~ NPOC_mg_C_per_L),
TN_mg_N_per_L = case_when(TN_mg_N_per_L == -9999 ~ NA,
str_detect(TN_mg_N_per_L, "LOD") ~ NA,
str_detect(TN_mg_N_per_L, "Below_0.1_ppm") ~
str_extract_part(TN_mg_N_per_L, "_ppm_Below_0.1_ppm", before = TRUE),
TRUE ~ TN_mg_N_per_L)) %>%
mutate(NPOC_mg_C_per_L = case_when(str_detect(NPOC_mg_C_per_L, "NPOC_") ~
str_extract(NPOC_mg_C_per_L, "(\\d+\\.\\d+)"),
NPOC_mg_C_per_L < 0 ~ NA,
TRUE ~ NPOC_mg_C_per_L),
TN_mg_N_per_L = case_when(str_detect(TN_mg_N_per_L, "TN_") ~
str_extract(TN_mg_N_per_L, "(\\d+\\.\\d+)"),
TN_mg_N_per_L < 0 ~ NA,
TRUE ~ TN_mg_N_per_L)) %>%
mutate(NPOC_mg_C_per_L = as.numeric(NPOC_mg_C_per_L),
TN_mg_N_per_L = as.numeric(TN_mg_N_per_L)) %>%
separate(col=Sample_ID, into = c("Sample_ID", "FilterSize"), sep="-")
SolidCN_clean <- essdive_template_cleaning(SolidCN) %>%
dplyr::rename(C_weight_per = `01463_C_percent`,
N_weight_per = `01472_N_percent`,
Parent_ID = Sample_Name) %>%
select(Parent_ID, C_weight_per, N_weight_per) %>%
filter(grepl('solid', Parent_ID)) %>%
mutate(Parent_ID = stringr::str_remove(Parent_ID ,"-solid"),
C_weight_per = case_when(C_weight_per == -9999 ~ NA,
str_detect(C_weight_per, "LOD") ~ NA,
TRUE ~ C_weight_per),
N_weight_per = case_when(N_weight_per == -9999 ~ NA,
str_detect(N_weight_per, "LOD") ~ NA,
TRUE ~ N_weight_per)) %>%
mutate(C_weight_per = as.numeric(C_weight_per),
N_weight_per = as.numeric(N_weight_per))
#negative values were replaced with NA, below LOD were replaced with NA, and above/below standard curve were replaced with the value
ICP_Solid_clean<- essdive_template_cleaning(ICP_Solid) %>%
dplyr::rename(Solid_Ca_mg_kg = '52718_Total_Ca_mg_per_kg',
Solid_Mg_mg_kg = '52719_Total_Mg_mg_per_kg',
Solid_P_mg_kg = 'Total_P_mg_per_kg',
Solid_K_mg_kg = '52720_Total_K_mg_per_kg',
Solid_Na_mg_kg = '52721_Total_Na_mg_per_kg',
Solid_S_mg_kg = '01466_Total_S_mg_per_kg',
Solid_Al_mg_kg = '01462_Total_Al_mg_per_kg',
Solid_Fe_mg_kg = '01464_Total_Fe_mg_per_kg',
Parent_ID = 'Sample_Name') %>%
mutate(Solid_Na_mg_kg = case_when(
Solid_Na_mg_kg == -9999 ~ NA,
str_detect(Solid_Na_mg_kg, "Negative_Value") ~ NA,
str_detect(Solid_Na_mg_kg, "Below_213.43_LOD") ~ NA,
str_detect(Solid_Na_mg_kg, "TotNa_Below") ~
str_extract(Solid_Na_mg_kg, "\\d+\\.?\\d*(?=_ppm_Final_Corrected)"),
str_detect(Solid_Na_mg_kg, "TotNa_Above") ~
str_extract(Solid_Na_mg_kg, "\\d+\\.?\\d*(?=_ppm_Final_Corrected)"),
TRUE ~ Solid_Na_mg_kg)) %>%
mutate_at(c(4:11), as.numeric) %>%
separate(col=Parent_ID, into = c("Parent_ID", "FilterSize"), sep="-") %>%
filter(grepl('solid', FilterSize)) %>%
select(Parent_ID, Solid_Ca_mg_kg, Solid_Mg_mg_kg, Solid_P_mg_kg, Solid_K_mg_kg, Solid_Na_mg_kg, Solid_S_mg_kg, Solid_Al_mg_kg, Solid_Fe_mg_kg)## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `Methods_Deviation =
## .Primitive("as.double")(Methods_Deviation)`.
## Caused by warning:
## ! NAs introduced by coercion#negative values were replaced with NA, below LOD were replaced with NA, and above/below standard curve were replaced with the value
ICP_Leachate_clean <- essdive_template_cleaning(ICP_Leachate) %>%
dplyr::rename(Leachate_Ca_mg_L = 'Total_Ca_mg_per_L',
Leachate_Mg_mg_L = 'Total_Mg_mg_per_L',
Leachate_P_mg_L = 'Total_P_mg_per_L',
Leachate_K_mg_L = 'Total_K_mg_per_L',
Leachate_Na_mg_L = 'Total_Na_mg_per_L',
Leachate_S_mg_L = 'Total_S_mg_per_L',
Leachate_Al_mg_L = 'Total_Al_mg_per_L',
Leachate_Fe_mg_L = 'Total_Fe_mg_per_L',
Sample_ID = 'Sample_Name') %>%
mutate(Leachate_Na_mg_L = case_when(
Leachate_Na_mg_L == -9999 ~ NA,
str_detect(Leachate_Na_mg_L, "Negative") ~ NA,
str_detect(Leachate_Na_mg_L, "LOD") ~ NA,
str_detect(Leachate_Na_mg_L, "TotNa_Below") ~
str_extract(Leachate_Na_mg_L, "\\d+\\.?\\d*(?=_ppm_Final_Corrected)"),
str_detect(Leachate_Na_mg_L, "TotNa_Above") ~
str_extract(Leachate_Na_mg_L, "\\d+\\.?\\d*(?=_ppm_Final_Corrected)"),
TRUE ~ Leachate_Na_mg_L)) %>%
mutate(Leachate_Al_mg_L = case_when(
Leachate_Na_mg_L == -9999 ~ NA,
str_detect(Leachate_Al_mg_L, "Negative") ~ NA,
str_detect(Leachate_Al_mg_L, "LOD") ~ NA,
str_detect(Leachate_Al_mg_L, "TotAl_Below") ~
str_extract(Leachate_Al_mg_L, "\\d+\\.?\\d*(?=_ppm_Final_Corrected)"),
str_detect(Leachate_Al_mg_L, "TotAl_Above") ~
str_extract(Leachate_Al_mg_L, "\\d+\\.?\\d*(?=_ppm_Final_Corrected)"),
TRUE ~ Leachate_Al_mg_L)) %>%
mutate(Leachate_Fe_mg_L = case_when(
Leachate_Fe_mg_L == -9999 ~ NA,
str_detect(Leachate_Fe_mg_L, "Negative") ~ NA,
str_detect(Leachate_Fe_mg_L, "LOD") ~ NA,
str_detect(Leachate_Fe_mg_L, "TotAl_Below") ~
str_extract(Leachate_Fe_mg_L, "\\d+\\.?\\d*(?=_ppm_Final_Corrected)"),
str_detect(Leachate_Fe_mg_L, "TotAl_Above") ~
str_extract(Leachate_Fe_mg_L, "\\d+\\.?\\d*(?=_ppm_Final_Corrected)"),
TRUE ~ Leachate_Fe_mg_L)) %>%
mutate_at(c(3:10), as.numeric) %>%
select(Sample_ID, Leachate_Ca_mg_L, Leachate_Mg_mg_L, Leachate_P_mg_L, Leachate_K_mg_L, Leachate_Na_mg_L, Leachate_S_mg_L, Leachate_Al_mg_L, Leachate_Fe_mg_L) %>%
filter(grepl('filt', Sample_ID)) %>%
separate(col=Sample_ID, into = c("Sample_ID", "FilterSize"), sep="-")
ICP_Leachate_Filt_clean<- ICP_Leachate_clean %>%
filter(str_detect(FilterSize, "filt0.")) %>%
dplyr::rename(Leachate_Ca_mg_L_filt0.7 = 'Leachate_Ca_mg_L',
Leachate_Mg_mg_L_filt0.7 = 'Leachate_Mg_mg_L',
Leachate_P_mg_L_filt0.7 = 'Leachate_P_mg_L',
Leachate_K_mg_L_filt0.7 = 'Leachate_K_mg_L',
Leachate_Na_mg_L_filt0.7 = 'Leachate_Na_mg_L',
Leachate_S_mg_L_filt0.7 = 'Leachate_S_mg_L',
Leachate_Al_mg_L_filt0.7 = 'Leachate_Al_mg_L',
Leachate_Fe_mg_L_filt0.7 = 'Leachate_Fe_mg_L') %>%
select(Sample_ID, Leachate_Ca_mg_L_filt0.7, Leachate_Mg_mg_L_filt0.7, Leachate_P_mg_L_filt0.7, Leachate_K_mg_L_filt0.7, Leachate_Na_mg_L_filt0.7, Leachate_S_mg_L_filt0.7, Leachate_Al_mg_L_filt0.7, Leachate_Fe_mg_L_filt0.7)
ICP_Leachate_Unfilt_clean<- ICP_Leachate_clean %>%
filter(str_detect(FilterSize, "unfilt")) %>%
dplyr::rename(Leachate_Ca_mg_L_unfilt = 'Leachate_Ca_mg_L',
Leachate_Mg_mg_L_unfilt = 'Leachate_Mg_mg_L',
Leachate_P_mg_L_unfilt = 'Leachate_P_mg_L',
Leachate_K_mg_L_unfilt = 'Leachate_K_mg_L',
Leachate_Na_mg_L_unfilt = 'Leachate_Na_mg_L',
Leachate_S_mg_L_unfilt = 'Leachate_S_mg_L',
Leachate_Al_mg_L_unfilt = 'Leachate_Al_mg_L',
Leachate_Fe_mg_L_unfilt = 'Leachate_Fe_mg_L') %>%
select(Sample_ID, Leachate_Ca_mg_L_unfilt, Leachate_Mg_mg_L_unfilt, Leachate_P_mg_L_unfilt, Leachate_K_mg_L_unfilt, Leachate_Na_mg_L_unfilt, Leachate_S_mg_L_unfilt, Leachate_Al_mg_L_unfilt, Leachate_Fe_mg_L_unfilt)
#negative values were replaced with zero, below LOD were replaced with NA, and above/below standard curve were replaced with the value (value is corrected for dilution and subtracted by unreactive P to correct for color interference)
#only used concentrations above the high standard if within 15% of it
MBD_Leachate_clean <- essdive_template_cleaning(MBD) %>%
dplyr::rename(Sample_ID = 'Sample_Name',
MBD_P_mg_per_L_filt0.7 = 'MBD_P_mg_per_L') %>%
mutate(MBD_P_mg_per_L_filt0.7 = case_when(
MBD_P_mg_per_L_filt0.7 == -9999 ~ NA,
str_detect(MBD_P_mg_per_L_filt0.7, "negative") ~ NA,
str_detect(MBD_P_mg_per_L_filt0.7, "LOD") ~ NA,
str_detect(Sample_ID, "50A") ~ NA, #hihger than 15% of high standard
str_detect(MBD_P_mg_per_L_filt0.7, "MBD_Below") ~
str_extract(MBD_P_mg_per_L_filt0.7, "\\d+\\.?\\d*(?=_ppm_Final_Corrected)"),
str_detect(MBD_P_mg_per_L_filt0.7, "MBD_Above") ~
str_extract(MBD_P_mg_per_L_filt0.7, "\\d+\\.?\\d*(?=_ppm_Final_Corrected)"),
TRUE ~ MBD_P_mg_per_L_filt0.7)) %>%
mutate_at(c(2), as.numeric) %>%
filter(grepl('filt', Sample_ID)) %>%
separate(col=Sample_ID, into = c("Sample_ID", "FilterSize"), sep="-") %>%
select(Sample_ID, MBD_P_mg_per_L_filt0.7) %>%
mutate_at(c(2), as.numeric) %>%
filter(Sample_ID != "BSLE_0073C") #remove 73C because below limit of quantification## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `Material = .Primitive("as.double")(Material)`.
## Caused by warning:
## ! NAs introduced by coercionICP_PNMR_clean <- essdive_template_cleaning(ICP_PNMR) %>%
dplyr::rename(Parent_ID = 'Sample_Name') %>%
mutate_at(c(2:8), as.numeric) %>%
filter(grepl('solid', Parent_ID)) %>%
separate(col=Parent_ID, into = c("Parent_ID", "FilterSize"), sep="-") %>%
select(Parent_ID, Total_Extractable_Ca_mg_per_kg, Total_Extractable_Mg_mg_per_kg, Total_Extractable_P_mg_per_kg, Total_Extractable_K_mg_per_kg, Total_Extractable_S_mg_per_kg, Total_Extractable_Al_mg_per_kg, Total_Extractable_Fe_mg_per_kg) %>%
mutate_at(c(2:8), as.numeric) ## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `Material = .Primitive("as.double")(Material)`.
## Caused by warning:
## ! NAs introduced by coercionXANES_clean <- XANES %>%
select(-where(~all(is.na(.)))) %>%
mutate(XANES_Po = rowSums(select(., contains("Po")), na.rm = TRUE)) %>%
mutate(XANES_Pi = rowSums(select(., contains("Pi")), na.rm = TRUE)) %>%
mutate(XANES_Fe = rowSums(select(., contains("Fe")), na.rm = TRUE)) %>%
mutate(XANES_Na = rowSums(select(., contains("Pi_Na") | contains("Po_Na")), na.rm = TRUE)) %>%
mutate(XANES_Ca = rowSums(select(., contains("Ca") | contains("Apatite")), na.rm = TRUE)) %>%
mutate(XANES_Al = rowSums(select(., contains("Al")), na.rm = TRUE)) %>%
mutate(XANES_Mg = rowSums(select(., contains("Mg") | contains("NH4_Mg")), na.rm = TRUE)) %>%
mutate(XANES_K = rowSums(select(., contains("K")), na.rm = TRUE)) %>%
mutate(XANES_Pi_Fe = rowSums(select(., contains("Pi_Fe")), na.rm = TRUE)) %>%
mutate(XANES_Pi_Na = rowSums(select(., contains("Pi_Na")), na.rm = TRUE)) %>%
mutate(XANES_Pi_Ca = rowSums(select(., contains("Pi_Ca") | contains("Apatite")), na.rm = TRUE)) %>%
mutate(XANES_Pi_Al = rowSums(select(., contains("Pi_Al")), na.rm = TRUE)) %>%
mutate(XANES_Pi_Mg = rowSums(select(., contains("Pi_Mg") | contains("NH4_Mg")), na.rm = TRUE)) %>%
mutate(XANES_Pi_K = rowSums(select(., contains("Pi_K")), na.rm = TRUE)) %>%
mutate(XANES_Po_Fe = rowSums(select(., contains("Po_Fe")), na.rm = TRUE)) %>%
mutate(XANES_Po_Na = rowSums(select(., contains("Po_Na")), na.rm = TRUE)) %>%
mutate(XANES_Po_Ca = rowSums(select(., contains("Po_Ca")), na.rm = TRUE)) %>%
mutate(XANES_Po_Al = rowSums(select(., contains("Po_Al")), na.rm = TRUE)) %>%
mutate(XANES_Po_Mg = rowSums(select(., contains("Po_Mg")), na.rm = TRUE)) %>%
mutate(XANES_Po_K = rowSums(select(., contains("Po_K")), na.rm = TRUE)) %>%
separate(col=Sample_Name, into = c("Parent_ID", "FilterSize"), sep="-") %>%
filter(str_detect(FilterSize, "solid"))
NMR_clean <- NMR %>%
separate(col=Sample_Name, into = c("Parent_ID", "FilterSize"), sep="-") data <- Meta %>%
full_join(SolidCN_clean, by = 'Parent_ID') %>%
full_join(ICP_Solid_clean, by = 'Parent_ID') %>%
full_join(ICP_PNMR_clean, by = 'Parent_ID') %>%
full_join(npoctdn_clean, by = 'Sample_ID') %>%
full_join(ICP_Leachate_Filt_clean, by = 'Sample_ID') %>%
full_join(ICP_Leachate_Unfilt_clean, by = 'Sample_ID') %>%
full_join(MBD_Leachate_clean, by = 'Sample_ID') %>%
full_join(pH_clean, by = 'Sample_ID') %>%
mutate_at(c(10:48), as.numeric) %>%
filter(!is.na(Solid_P_mg_kg))## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `FilterSize = .Primitive("as.double")(FilterSize)`.
## Caused by warning:
## ! NAs introduced by coerciondata$Burn_Severity=factor(data$Burn_Severity, levels= c("Unburned", "Low", "Moderate", "High"))
#create a dataframe for only solids figures and stats to prevent having reps count as extra samples; add NMR and XANES data
data_solids <- data %>%
filter(str_detect(Sample_ID, "A")) %>%
select(1:28) %>%
left_join(NMR_clean) %>%
left_join(XANES_clean)## Joining with `by = join_by(Parent_ID)`
## Joining with `by = join_by(Parent_ID, FilterSize)`data_solids$Burn_Severity=factor(data_solids$Burn_Severity, levels= c("Unburned", "Low", "Moderate", "High"))
#keep XANES separate too
XANES_clean<- Meta %>%
full_join(XANES_clean, by = 'Parent_ID') %>%
filter(str_detect(Sample_ID, "A")) %>%
mutate(Burn_Treatment = as.character(Burn_Treatment)) %>%
filter(Burn_Treatment %in% c("burn table" , "raw")) %>%
filter(Land_Coverage_Category %in% c("Douglas-fir forest" , "Sagebrush shrubland")) %>%
filter(!is.na(FilterSize))####Additional Calcs
#NMR extraction efficiency
data_solids<-data_solids %>%
mutate(NMR_EE_P_perc=Total_Extractable_P_mg_per_kg/Solid_P_mg_kg*100)
#convert C and N of solids from percent to mg/kg
data_solids<-data_solids %>%
mutate(Solid_C_mg_kg = (C_weight_per*10000),
Solid_N_mg_kg = (N_weight_per*10000))
#convert to P relative to solid mass and water volume for Molybdate
data<-data %>%
mutate(MBD_P_mg_kgchar = MBD_P_mg_per_L_filt0.7*Leachate_vol_mL/Char_leached_g*1000/1000)
#normalize by P concentration in solid char for Molybdate
data<-data %>%
mutate(MBD_P_percentofsolid = MBD_P_mg_kgchar/Solid_P_mg_kg*100)
#calculate particulate concentrations as the difference between filtered and unfiltered
data<-data %>%
mutate(Leachate_Ca_mg_L_part = Leachate_Ca_mg_L_unfilt - Leachate_Ca_mg_L_filt0.7) %>%
mutate(Leachate_Mg_mg_L_part = Leachate_Mg_mg_L_unfilt - Leachate_Mg_mg_L_filt0.7) %>%
mutate(Leachate_P_mg_L_part = Leachate_P_mg_L_unfilt - Leachate_P_mg_L_filt0.7) %>%
mutate(Leachate_K_mg_L_part = Leachate_K_mg_L_unfilt - Leachate_K_mg_L_filt0.7) %>%
mutate(Leachate_Na_mg_L_part = Leachate_Na_mg_L_unfilt - Leachate_Na_mg_L_filt0.7) %>%
mutate(Leachate_S_mg_L_part = Leachate_S_mg_L_unfilt - Leachate_S_mg_L_filt0.7) %>%
mutate(Leachate_Al_mg_L_part = Leachate_Al_mg_L_unfilt - Leachate_Al_mg_L_filt0.7) %>%
mutate(Leachate_Fe_mg_L_part = Leachate_Fe_mg_L_unfilt - Leachate_Fe_mg_L_filt0.7)
#convert to P relative to solid mass and water volume for the soluble and particulate fractions
data<-data %>%
mutate(Leachate_Ca_mg_kgchar_part = Leachate_Ca_mg_L_part*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
mutate(Leachate_Mg_mg_kgchar_part = Leachate_Mg_mg_L_part*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
mutate(Leachate_P_mg_kgchar_part = Leachate_P_mg_L_part*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
mutate(Leachate_K_mg_kgchar_part = Leachate_K_mg_L_part*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
mutate(Leachate_Na_mg_kgchar_part = Leachate_Na_mg_L_part*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
mutate(Leachate_S_mg_kgchar_part = Leachate_S_mg_L_part*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
mutate(Leachate_Al_mg_kgchar_part = Leachate_Al_mg_L_part*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
mutate(Leachate_Fe_mg_kgchar_part = Leachate_Fe_mg_L_part*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
mutate(Leachate_Ca_mg_kgchar_filt0.7 = Leachate_Ca_mg_L_filt0.7*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
mutate(Leachate_Mg_mg_kgchar_filt0.7 = Leachate_Mg_mg_L_filt0.7*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
mutate(Leachate_P_mg_kgchar_filt0.7 = Leachate_P_mg_L_filt0.7*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
mutate(Leachate_K_mg_kgchar_filt0.7 = Leachate_K_mg_L_filt0.7*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
mutate(Leachate_Na_mg_kgchar_filt0.7 = Leachate_Na_mg_L_filt0.7*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
mutate(Leachate_S_mg_kgchar_filt0.7 = Leachate_S_mg_L_filt0.7*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
mutate(Leachate_Al_mg_kgchar_filt0.7 = Leachate_Al_mg_L_filt0.7*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
mutate(Leachate_Fe_mg_kgchar_filt0.7 = Leachate_Fe_mg_L_filt0.7*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
mutate(Leachate_Ca_mg_kgchar_unfilt = Leachate_Ca_mg_L_unfilt*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
mutate(Leachate_Mg_mg_kgchar_unfilt = Leachate_Mg_mg_L_unfilt*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
mutate(Leachate_P_mg_kgchar_unfilt = Leachate_P_mg_L_unfilt*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
mutate(Leachate_K_mg_kgchar_unfilt = Leachate_K_mg_L_unfilt*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
mutate(Leachate_Na_mg_kgchar_unfilt = Leachate_Na_mg_L_unfilt*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
mutate(Leachate_S_mg_kgchar_unfilt = Leachate_S_mg_L_unfilt*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
mutate(Leachate_Al_mg_kgchar_unfilt = Leachate_Al_mg_L_unfilt*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
mutate(Leachate_Fe_mg_kgchar_unfilt = Leachate_Fe_mg_L_unfilt*Leachate_vol_mL/Char_leached_g*1000/1000)
#convert solid elemental concentrations from to g/kg
data_solids<-data_solids %>%
mutate(Solid_P_g_kg = Solid_P_mg_kg/1000) %>%
mutate(Solid_Ca_g_kg = Solid_Ca_mg_kg/1000) %>%
mutate(Solid_Fe_g_kg = Solid_Fe_mg_kg/1000) %>%
mutate(Solid_Al_g_kg = Solid_Al_mg_kg/1000) %>%
mutate(Solid_K_g_kg = Solid_K_mg_kg/1000) %>%
mutate(Solid_Mg_g_kg = Solid_Mg_mg_kg/1000) %>%
mutate(Solid_Na_g_kg = Solid_Na_mg_kg/1000) %>%
mutate(Solid_S_g_kg = Solid_S_mg_kg/1000) %>%
mutate(Solid_C_g_kg = C_weight_per/100*1000) %>%
mutate(Solid_N_g_kg = N_weight_per/100*1000)
#normalize by P concentration in solid char fix
data<-data %>%
mutate(Leachate_Ca_perc_solid_part = Leachate_Ca_mg_kgchar_part/Solid_P_mg_kg*100) %>%
mutate(Leachate_Mg_perc_solid_part = Leachate_Mg_mg_kgchar_part/Solid_P_mg_kg*100) %>%
mutate(Leachate_P_perc_solid_part = Leachate_P_mg_kgchar_part/Solid_P_mg_kg*100) %>%
mutate(Leachate_K_perc_solid_part = Leachate_K_mg_kgchar_part/Solid_P_mg_kg*100) %>%
mutate(Leachate_Na_perc_solid_part = Leachate_Na_mg_kgchar_part/Solid_P_mg_kg*100) %>%
mutate(Leachate_S_perc_solid_part = Leachate_S_mg_kgchar_part/Solid_P_mg_kg*100) %>%
mutate(Leachate_Al_perc_solid_part = Leachate_Al_mg_kgchar_part/Solid_P_mg_kg*100) %>%
mutate(Leachate_Fe_perc_solid_part = Leachate_Fe_mg_kgchar_part/Solid_P_mg_kg*100) %>%
mutate(Leachate_Ca_perc_solid_filt0.7 = Leachate_Ca_mg_kgchar_filt0.7/Solid_P_mg_kg*100) %>%
mutate(Leachate_Mg_perc_solid_filt0.7 = Leachate_Mg_mg_kgchar_filt0.7/Solid_P_mg_kg*100) %>%
mutate(Leachate_P_perc_solid_filt0.7 = Leachate_P_mg_kgchar_filt0.7/Solid_P_mg_kg*100) %>%
mutate(Leachate_K_perc_solid_filt0.7 = Leachate_K_mg_kgchar_filt0.7/Solid_P_mg_kg*100) %>%
mutate(Leachate_Na_perc_solid_filt0.7 = Leachate_Na_mg_kgchar_filt0.7/Solid_P_mg_kg*100) %>%
mutate(Leachate_S_perc_solid_filt0.7 = Leachate_S_mg_kgchar_filt0.7/Solid_P_mg_kg*100) %>%
mutate(Leachate_Al_perc_solid_filt0.7 = Leachate_Al_mg_kgchar_filt0.7/Solid_P_mg_kg*100) %>%
mutate(Leachate_Fe_perc_solid_filt0.7 = Leachate_Fe_mg_kgchar_filt0.7/Solid_P_mg_kg*100)
#molar stoichiometry of P:element in solid phase fix
test<-data_solids %>%
mutate(Fe.P.MolarRatio=(Solid_Fe_mg_kg/55.845/1000)/(Solid_P_mg_kg/30.97/1000))%>%
mutate(K.P.MolarRatio=(Solid_K_mg_kg/39.0983/1000)/(Solid_P_mg_kg/30.97/1000))%>%
mutate(Na.P.MolarRatio=(Solid_Na_mg_kg/22.989/1000)/(Solid_P_mg_kg/30.97/1000))%>%
mutate(Al.P.MolarRatio=(Solid_Al_mg_kg/26.982/1000)/(Solid_P_mg_kg/30.97/1000))%>%
mutate(Mg.P.MolarRatio=(Solid_Mg_mg_kg/24.305/1000)/(Solid_P_mg_kg/30.97/1000))%>%
mutate(Ca.P.MolarRatio=(Solid_Ca_mg_kg/40.078/1000)/(Solid_P_mg_kg/30.97/1000))%>%
mutate(S.P.MolarRatio=(Solid_S_mg_kg/32.065/1000)/(Solid_P_mg_kg/30.97/1000))
#molar stoichiometry of P:element in particulate phase
#molar stoichiometry of P:element in <0.7 phase
#save csv of final collated data with manipulations
write.csv(data, "../data/data.csv", row.names = FALSE) #all data
write.csv(data_solids, "../data/data_solids.csv", row.names = FALSE) #solids only##Figures and Stats
#import data so you don't have to run above code
data <- read_csv("../data/data.csv")## Rows: 57 Columns: 99
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (7): Parent_ID, Sample_ID, Land_Coverage_Category, Burn_Treatment, Burn...
## dbl (91): Char_leached_g, Leachate_vol_mL, Burn_Duration, Char_Max_Temp, C_w...
## lgl (1): FilterSize
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.data_solids <- read_csv("../data/data_solids.csv")## Rows: 19 Columns: 120
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (10): Parent_ID, Sample_ID, Land_Coverage_Category, Burn_Treatment, Bur...
## dbl (110): Char_leached_g, Leachate_vol_mL, Burn_Duration, Char_Max_Temp, C_...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.####summary stats
#solid elemental concentrations, NMR, and XANES summary
summary_solids_conc <- data_solids %>%
dplyr::group_by(Burn_Severity, Land_Coverage_Category) %>%
dplyr::summarize(
P_Mean = mean(Solid_P_g_kg, na.rm = TRUE),
P_Standard_Deviation = sd(Solid_P_g_kg, na.rm = TRUE),
Ca_Mean = mean(Solid_Ca_g_kg, na.rm = TRUE),
Ca_Standard_Deviation = sd(Solid_Ca_g_kg, na.rm = TRUE),
Fe_Mean = mean(Solid_Fe_g_kg, na.rm = TRUE),
Fe_Standard_Deviation = sd(Solid_Fe_g_kg, na.rm = TRUE),
Al_Mean = mean(Solid_Al_g_kg, na.rm = TRUE),
Al_Standard_Deviation = sd(Solid_Al_g_kg, na.rm = TRUE),
K_Mean = mean(Solid_K_g_kg, na.rm = TRUE),
K_Standard_Deviation = sd(Solid_K_g_kg, na.rm = TRUE),
Mg_Mean = mean(Solid_Mg_g_kg, na.rm = TRUE),
Mg_Standard_Deviation = sd(Solid_Mg_g_kg, na.rm = TRUE),
Na_Mean = mean(Solid_Na_g_kg, na.rm = TRUE),
Na_Standard_Deviation = sd(Solid_Na_g_kg, na.rm = TRUE),
S_Mean = mean(Solid_S_g_kg, na.rm = TRUE),
S_Standard_Deviation = sd(Solid_S_g_kg, na.rm = TRUE),
C_Mean = mean(Solid_C_g_kg, na.rm = TRUE),
C_Standard_Deviation = sd(Solid_C_g_kg, na.rm = TRUE),
N_Mean = mean(Solid_N_g_kg, na.rm = TRUE),
N_Standard_Deviation = sd(Solid_N_g_kg, na.rm = TRUE),
Ortho_Mean = mean(Ortho, na.rm = TRUE),
Ortho_Standard_Deviation = sd(Ortho, na.rm = TRUE),
Mono_Mean = mean(Mono, na.rm = TRUE),
Mono_Standard_Deviation = sd(Mono, na.rm = TRUE),
Di_Mean = mean(Di, na.rm = TRUE),
Di_Standard_Deviation = sd(Di, na.rm = TRUE),
Pyro_Mean = mean(Pyro, na.rm = TRUE),
Pyro_Standard_Deviation = sd(Pyro, na.rm = TRUE),
XANES_Po_Mean = mean(XANES_Po, na.rm = TRUE),
XANES_Po_Standard_Deviation = sd(XANES_Po, na.rm = TRUE),
XANES_Pi_Mean = mean(XANES_Pi, na.rm = TRUE),
XANES_Pi_Standard_Deviation = sd(XANES_Pi, na.rm = TRUE),
XANES_Fe_Mean = mean(XANES_Fe, na.rm = TRUE),
XANES_Fe_Standard_Deviation = sd(XANES_Fe, na.rm = TRUE),
XANES_Na_Mean = mean(XANES_Na, na.rm = TRUE),
XANES_Na_Standard_Deviation = sd(XANES_Na, na.rm = TRUE),
XANES_Ca_Mean = mean(XANES_Ca, na.rm = TRUE),
XANES_Ca_Standard_Deviation = sd(XANES_Ca, na.rm = TRUE),
XANES_Al_Mean = mean(XANES_Al, na.rm = TRUE),
XANES_Al_Standard_Deviation = sd(XANES_Al, na.rm = TRUE),
XANES_Mg_Mean = mean(XANES_Mg, na.rm = TRUE),
XANES_Mg_Standard_Deviation = sd(XANES_Mg, na.rm = TRUE),
XANES_K_Mean = mean(XANES_K, na.rm = TRUE),
XANES_K_Standard_Deviation = sd(XANES_K, na.rm = TRUE),
XANES_Pi_Fe_Mean = mean(XANES_Pi_Fe, na.rm = TRUE),
XANES_Pi_Fe_Standard_Deviation = sd(XANES_Pi_Fe, na.rm = TRUE),
XANES_Pi_Na_Mean = mean(XANES_Pi_Na, na.rm = TRUE),
XANES_Pi_Na_Standard_Deviation = sd(XANES_Pi_Na, na.rm = TRUE),
XANES_Pi_Ca_Mean = mean(XANES_Pi_Ca, na.rm = TRUE),
XANES_Pi_Ca_Standard_Deviation = sd(XANES_Pi_Ca, na.rm = TRUE),
XANES_Pi_Al_Mean = mean(XANES_Pi_Al, na.rm = TRUE),
XANES_Pi_Al_Standard_Deviation = sd(XANES_Pi_Al, na.rm = TRUE),
XANES_Pi_Mg_Mean = mean(XANES_Pi_Mg, na.rm = TRUE),
XANES_Pi_Mg_Standard_Deviation = sd(XANES_Pi_Mg, na.rm = TRUE),
XANES_Pi_K_Mean = mean(XANES_Pi_K, na.rm = TRUE),
XANES_Pi_K_Standard_Deviation = sd(XANES_Pi_K, na.rm = TRUE),
XANES_Po_Fe_Mean = mean(XANES_Po_Fe, na.rm = TRUE),
XANES_Po_Fe_Standard_Deviation = sd(XANES_Po_Fe, na.rm = TRUE),
XANES_Po_Na_Mean = mean(XANES_Po_Na, na.rm = TRUE),
XANES_Po_Na_Standard_Deviation = sd(XANES_Po_Na, na.rm = TRUE),
XANES_Po_Ca_Mean = mean(XANES_Po_Ca, na.rm = TRUE),
XANES_Po_Ca_Standard_Deviation = sd(XANES_Po_Ca, na.rm = TRUE),
XANES_Po_Al_Mean = mean(XANES_Po_Al, na.rm = TRUE),
XANES_Po_Al_Standard_Deviation = sd(XANES_Po_Al, na.rm = TRUE),
XANES_Po_Mg_Mean = mean(XANES_Po_Mg, na.rm = TRUE),
XANES_Po_Mg_Standard_Deviation = sd(XANES_Po_Mg, na.rm = TRUE),
XANES_Po_K_Mean = mean(XANES_Po_K, na.rm = TRUE),
XANES_Po_K_Standard_Deviation = sd(XANES_Po_K, na.rm = TRUE),
Count = n()
)## `summarise()` has grouped output by 'Burn_Severity'. You can override using the
## `.groups` argument.write.csv(summary_solids_conc, "../data/summary_solids_conc.csv") #save this for a table
#NMR extraction efficiency
summary_solids <- data_solids %>%
dplyr::group_by(Burn_Severity, Land_Coverage_Category) %>%
dplyr::summarize(
Mean = mean(NMR_EE_P_perc, na.rm = TRUE),
Standard_Deviation = sd(NMR_EE_P_perc, na.rm = TRUE),
Count = n()
)## `summarise()` has grouped output by 'Burn_Severity'. You can override using the
## `.groups` argument.#leachate concentrations
P.Conc.leach<- data %>%
dplyr::group_by(Burn_Severity, Land_Coverage_Category) %>%
dplyr::summarize(
Leachate_P_mg_kgchar_part_Mean = mean(Leachate_P_mg_kgchar_part, na.rm = TRUE),
Leachate_P_mg_kgchar_part_Standard_Deviation = sd(Leachate_P_mg_kgchar_part, na.rm = TRUE),
Count = n(),
Leachate_P_mg_kgchar_filt0.7_Mean = mean(Leachate_P_mg_kgchar_filt0.7, na.rm = TRUE),
Leachate_P_mg_kgchar_filt0.7_Standard_Deviation = sd(Leachate_P_mg_kgchar_filt0.7, na.rm = TRUE),
Leachate_P_mg_kgchar_unfilt_Mean = mean(Leachate_P_mg_kgchar_unfilt, na.rm = TRUE),
Leachate_P_mg_kgchar_unfilt_Standard_Deviation = sd(Leachate_P_mg_kgchar_unfilt, na.rm = TRUE)
)## `summarise()` has grouped output by 'Burn_Severity'. You can override using the
## `.groups` argument.####Molybdate
#MBD mg/L figure
ggplot(data, aes(x=Burn_Severity, y=MBD_P_mg_per_L_filt0.7, fill=Land_Coverage_Category))+geom_boxplot(show.legend=FALSE)+theme_classic()+ylab(expression(paste(Soluble~Reactive~P~(mg~P~L^-1), sep="")))+ xlab("Burn Severity")+scale_fill_manual(values=c("#346648", "#dbb40d"))+facet_grid(.~Land_Coverage_Category)+theme(text = element_text(size = 15))## Warning: Removed 3 rows containing non-finite values (`stat_boxplot()`).
#MBD mg/kg figure; likely use this one for publication
MBD.BS.fig<-ggplot(data, aes(x=Burn_Severity, y=MBD_P_mg_kgchar/1000, fill=Land_Coverage_Category))+geom_boxplot(show.legend=FALSE)+theme_classic()+ylab(expression(paste(Soluble~Reactive~P~(g~P~kg^-1), sep="")))+ xlab("Burn Severity")+scale_fill_manual(values=c("#346648", "#dbb40d"))+facet_grid(.~Land_Coverage_Category)+theme(text = element_text(size = 15))
cowplot::save_plot("../figures/MBD.BS.fig.pdf", MBD.BS.fig, ncol = 1, nrow = 1, base_aspect_ratio= 2:1, dpi=72)## Warning: Removed 3 rows containing non-finite values (`stat_boxplot()`).#MBD as percent of solid P concentration figure
ggplot(data, aes(x=Burn_Severity, y=MBD_P_percentofsolid, fill=Land_Coverage_Category))+geom_boxplot(show.legend=FALSE)+theme_classic()+ylab("Soluble Reactive P (%)")+ xlab("Burn Severity")+scale_fill_manual(values=c("#346648", "#dbb40d"))+facet_grid(.~Land_Coverage_Category)+theme(text = element_text(size = 15))## Warning: Removed 3 rows containing non-finite values (`stat_boxplot()`).
#MBD g/kg stats; likely use this one for publication
model<-lmer(log10(MBD_P_mg_kgchar/1000)~Burn_Severity*Land_Coverage_Category+(1|Parent_ID), data=data) #if want to use maximum likelihood , REML=FALSE## fixed-effect model matrix is rank deficient so dropping 1 column / coefficientmodel.BS<-lmer(log10(MBD_P_mg_kgchar/1000)~Land_Coverage_Category+(1|Parent_ID), data=data)
model.LC<-lmer(log10(MBD_P_mg_kgchar/1000)~Burn_Severity+(1|Parent_ID), data=data)
model.int<-lmer(log10(MBD_P_mg_kgchar/1000)~Burn_Severity+Land_Coverage_Category+(1|Parent_ID), data=data)
anova(model,model.BS) #not sig## refitting model(s) with ML (instead of REML)## Data: data
## Models:
## model.BS: log10(MBD_P_mg_kgchar/1000) ~ Land_Coverage_Category + (1 | Parent_ID)
## model: log10(MBD_P_mg_kgchar/1000) ~ Burn_Severity * Land_Coverage_Category + (1 | Parent_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## model.BS 4 30.181 38.137 -11.0905 22.181
## model 9 32.325 50.226 -7.1625 14.325 7.856 5 0.1644anova(model,model.LC) #sig## refitting model(s) with ML (instead of REML)## Data: data
## Models:
## model.LC: log10(MBD_P_mg_kgchar/1000) ~ Burn_Severity + (1 | Parent_ID)
## model: log10(MBD_P_mg_kgchar/1000) ~ Burn_Severity * Land_Coverage_Category + (1 | Parent_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## model.LC 6 37.291 49.225 -12.6453 25.291
## model 9 32.325 50.226 -7.1625 14.325 10.966 3 0.01191 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(model,model.int) #sig## refitting model(s) with ML (instead of REML)## Data: data
## Models:
## model.int: log10(MBD_P_mg_kgchar/1000) ~ Burn_Severity + Land_Coverage_Category + (1 | Parent_ID)
## model: log10(MBD_P_mg_kgchar/1000) ~ Burn_Severity * Land_Coverage_Category + (1 | Parent_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## model.int 7 35.410 49.333 -10.7051 21.410
## model 9 32.325 50.226 -7.1625 14.325 7.0852 2 0.02894 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#check VIF
vif(model)## GVIF Df GVIF^(1/(2*Df))
## Burn_Severity 2.567759 3 1.170197
## Land_Coverage_Category 5.502341 1 2.345707
## Burn_Severity:Land_Coverage_Category 9.075768 2 1.735685#model assumptions
hist(log10(data$MBD_P_mg_kgchar))
#linearity
ggqqplot(residuals(model))
#check normality of residuals
shapiro.test(log10(data$MBD_P_mg_kgchar)) #p > 0.05 so can assume normality ##
## Shapiro-Wilk normality test
##
## data: log10(data$MBD_P_mg_kgchar)
## W = 0.95892, p-value = 0.06167#check homogeneity of variance assumption (homoscedastcity)
plot(fitted(model),residuals(model))
performance::check_heteroscedasticity(model) #p > 0.05 so can assume homoscedastic## OK: Error variance appears to be homoscedastic (p = 0.632).#post hoc tests for interaction term
emmeans::lsmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Burn_Severity="Unburned"),adjust="tukey")## $lsmeans
## Burn_Severity Land_Coverage_Category lsmean SE df lower.CL upper.CL
## Unburned Douglas-fir forest -1.183 0.221 11.7 -1.67 -0.6994
## Unburned Sagebrush shrubland -0.701 0.313 11.7 -1.39 -0.0163
##
## Degrees-of-freedom method: kenward-roger
## Results are given on the log10 (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE
## (Unburned Douglas-fir forest) - Unburned Sagebrush shrubland -0.483 0.384
## df t.ratio p.value
## 11.7 -1.258 0.2328
##
## Degrees-of-freedom method: kenward-roger
## Results are given on the log10 (not the response) scale.emmeans::lsmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Burn_Severity="Low"),adjust="tukey")## $lsmeans
## Burn_Severity Land_Coverage_Category lsmean SE df lower.CL upper.CL
## Low Douglas-fir forest -1.266 0.140 11.7 -1.57 -0.960
## Low Sagebrush shrubland -0.551 0.221 11.7 -1.03 -0.067
##
## Degrees-of-freedom method: kenward-roger
## Results are given on the log10 (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio
## (Low Douglas-fir forest) - Low Sagebrush shrubland -0.715 0.262 11.7 -2.729
## p.value
## 0.0187
##
## Degrees-of-freedom method: kenward-roger
## Results are given on the log10 (not the response) scale.emmeans::lsmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Burn_Severity="Moderate"),adjust="tukey")## $lsmeans
## Burn_Severity Land_Coverage_Category lsmean SE df lower.CL upper.CL
## Moderate Douglas-fir forest -1.05 0.181 11.7 -1.45 -0.655
## Moderate Sagebrush shrubland -1.23 0.221 11.7 -1.71 -0.742
##
## Degrees-of-freedom method: kenward-roger
## Results are given on the log10 (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE
## (Moderate Douglas-fir forest) - Moderate Sagebrush shrubland 0.176 0.286
## df t.ratio p.value
## 11.7 0.616 0.5500
##
## Degrees-of-freedom method: kenward-roger
## Results are given on the log10 (not the response) scale.emmeans::lsmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Land_Coverage_Category="Douglas-fir forest"),adjust="tukey")## $lsmeans
## Burn_Severity Land_Coverage_Category lsmean SE df lower.CL upper.CL
## High Douglas-fir forest -1.07 0.161 13.1 -1.42 -0.726
## Low Douglas-fir forest -1.27 0.140 11.7 -1.57 -0.960
## Moderate Douglas-fir forest -1.05 0.181 11.7 -1.45 -0.655
## Unburned Douglas-fir forest -1.18 0.221 11.7 -1.67 -0.699
##
## Degrees-of-freedom method: kenward-roger
## Results are given on the log10 (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE
## (High Douglas-fir forest) - (Low Douglas-fir forest) 0.1918 0.214
## (High Douglas-fir forest) - (Moderate Douglas-fir forest) -0.0238 0.242
## (High Douglas-fir forest) - (Unburned Douglas-fir forest) 0.1093 0.274
## (Low Douglas-fir forest) - (Moderate Douglas-fir forest) -0.2156 0.229
## (Low Douglas-fir forest) - (Unburned Douglas-fir forest) -0.0825 0.262
## (Moderate Douglas-fir forest) - (Unburned Douglas-fir forest) 0.1331 0.286
## df t.ratio p.value
## 12.4 0.898 0.8062
## 12.3 -0.098 0.9996
## 12.1 0.399 0.9775
## 11.7 -0.943 0.7830
## 11.7 -0.315 0.9886
## 11.7 0.466 0.9651
##
## Degrees-of-freedom method: kenward-roger
## Results are given on the log10 (not the response) scale.
## P value adjustment: tukey method for comparing a family of 4 estimatesemmeans::lsmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Land_Coverage_Category="Sagebrush shrubland"),adjust="tukey")## $lsmeans
## Burn_Severity Land_Coverage_Category lsmean SE df lower.CL upper.CL
## High Sagebrush shrubland nonEst NA NA NA NA
## Low Sagebrush shrubland -0.551 0.221 11.7 -1.03 -0.0670
## Moderate Sagebrush shrubland -1.226 0.221 11.7 -1.71 -0.7423
## Unburned Sagebrush shrubland -0.701 0.313 11.7 -1.39 -0.0163
##
## Degrees-of-freedom method: kenward-roger
## Results are given on the log10 (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE
## High Sagebrush shrubland - Low Sagebrush shrubland nonEst NA
## High Sagebrush shrubland - Moderate Sagebrush shrubland nonEst NA
## High Sagebrush shrubland - Unburned Sagebrush shrubland nonEst NA
## Low Sagebrush shrubland - Moderate Sagebrush shrubland 0.675 0.313
## Low Sagebrush shrubland - Unburned Sagebrush shrubland 0.150 0.384
## Moderate Sagebrush shrubland - Unburned Sagebrush shrubland -0.526 0.384
## df t.ratio p.value
## NA NA NA
## NA NA NA
## NA NA NA
## 11.7 2.157 0.1908
## 11.7 0.391 0.9788
## 11.7 -1.370 0.5396
##
## Degrees-of-freedom method: kenward-roger
## Results are given on the log10 (not the response) scale.
## P value adjustment: tukey method for comparing a family of 4 estimates####ICP Leachates
#convert to percent of total P in unfiltered sample
#ICP$P.filt.percentofunfilt<-(ICP$Total_P_mg_per_L_filt/ICP$Total_P_mg_per_L_unfilt)*100
#ICP$P.part.percentofunfilt<-(ICP$Total_P_mg_per_L_part/ICP$Total_P_mg_per_L_unfilt)*100
#total P in unfiltered leachate
Leachate.total.P.BS.boxplot <- ggplot(data, aes(x=Burn_Severity, y=Leachate_P_mg_kgchar_unfilt, fill=Land_Coverage_Category))+geom_boxplot(show.legend=FALSE)+theme_classic()+ylab(expression(paste(Total~P~(mg~P~kg^-1), sep="")))+ xlab("Burn Severity")+scale_fill_manual(values=c("#346648", "#dbb40d"))+facet_grid(.~Land_Coverage_Category)+theme(text = element_text(size = 15))
cowplot::save_plot("../figures/Leachate.total.P.BS.boxplot.pdf", Leachate.total.P.BS.boxplot, ncol = 1, nrow = 1, base_aspect_ratio= 2:1, dpi=72)
Leachate.filt0.7.total.P.BS.boxplot <- ggplot(data, aes(x=Burn_Severity, y=Leachate_P_mg_kgchar_filt0.7, fill=Land_Coverage_Category))+geom_boxplot(show.legend=FALSE)+theme_classic()+ylab(expression(paste(Total~Soluble~P~(mg~P~kg^-1), sep="")))+ xlab("Burn Severity")+scale_fill_manual(values=c("#346648", "#dbb40d"))+facet_grid(.~Land_Coverage_Category)+theme(text = element_text(size = 15))
cowplot::save_plot("../figures/Leachate.filt0.7.total.P.BS.boxplot.pdf", Leachate.filt0.7.total.P.BS.boxplot, ncol = 1, nrow = 1, base_aspect_ratio= 2:1, dpi=72)
Leachate.unfilt.total.P.BS.boxplot <- ggplot(data, aes(x=Burn_Severity, y=Leachate_P_mg_kgchar_part, fill=Land_Coverage_Category))+geom_boxplot(show.legend=FALSE)+theme_classic()+ylab(expression(paste(Total~Particulate~P~(mg~P~kg^-1), sep="")))+ xlab("Burn Severity")+scale_fill_manual(values=c("#346648", "#dbb40d"))+facet_grid(.~Land_Coverage_Category)+theme(text = element_text(size = 15))
cowplot::save_plot("../figures/Leachate.unfilt.total.P.BS.boxplot.pdf", Leachate.unfilt.total.P.BS.boxplot, ncol = 1, nrow = 1, base_aspect_ratio= 2:1, dpi=72)
#stacked bar chart
P.Conc.leach.bar <- P.Conc.leach %>%
mutate(sd.high=Leachate_P_mg_kgchar_unfilt_Mean+Leachate_P_mg_kgchar_unfilt_Standard_Deviation) %>%
mutate(sd.low=Leachate_P_mg_kgchar_unfilt_Mean-Leachate_P_mg_kgchar_unfilt_Standard_Deviation)
P.Conc.leach_melt<-melt(P.Conc.leach.bar, id=c("Burn_Severity", "Land_Coverage_Category","sd.low", "sd.high")) %>%
filter(variable %in% c("Leachate_P_mg_kgchar_part_Mean" , "Leachate_P_mg_kgchar_filt0.7_Mean"))
Leachate_stacked_bar <- ggplot(subset(P.Conc.leach_melt, variable=="Leachate_P_mg_kgchar_part_Mean" | variable=="Leachate_P_mg_kgchar_filt0.7_Mean"), aes(x = Burn_Severity, y = value, fill=variable)) + geom_bar(stat = "identity") + theme_classic()+ylab(expression(paste(Leachate~Mean~Total~P~(mg~P~kg^-1), sep="")))+ xlab("Burn Severity")+scale_fill_manual(values=c("#2274A5", "#E7DFC6"), labels = c("Particulate", "Soluble"), name = "")+facet_grid(.~Land_Coverage_Category)+theme(text = element_text(size = 15)) + geom_errorbar(aes(ymax=sd.high, ymin=sd.low), width = 0.2, size = 1, color = "black")## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.cowplot::save_plot("../figures/Leachate_stacked_bar.pdf", Leachate_stacked_bar, ncol = 1, nrow = 1, base_aspect_ratio= 2:1, dpi=72)
#ICP unfiltered total P g/kg stats; likely use this one for publication
model<-lmer(log10(Leachate_P_mg_kgchar_unfilt/1000)~Burn_Severity*Land_Coverage_Category+(1|Parent_ID), data=data)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficientmodel.BS<-lmer(log10(Leachate_P_mg_kgchar_unfilt/1000)~Land_Coverage_Category+(1|Parent_ID), data=data)
model.LC<-lmer(log10(Leachate_P_mg_kgchar_unfilt/1000)~Burn_Severity+(1|Parent_ID), data=data)
model.int<-lmer(log10(Leachate_P_mg_kgchar_unfilt/1000)~Burn_Severity+Land_Coverage_Category+(1|Parent_ID), data=data)
anova(model,model.BS)## refitting model(s) with ML (instead of REML)## Data: data
## Models:
## model.BS: log10(Leachate_P_mg_kgchar_unfilt/1000) ~ Land_Coverage_Category + (1 | Parent_ID)
## model: log10(Leachate_P_mg_kgchar_unfilt/1000) ~ Burn_Severity * Land_Coverage_Category + (1 | Parent_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## model.BS 4 44.894 53.066 -18.4470 36.894
## model 9 24.615 43.003 -3.3077 6.615 30.279 5 1.3e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(model,model.LC)## refitting model(s) with ML (instead of REML)## Data: data
## Models:
## model.LC: log10(Leachate_P_mg_kgchar_unfilt/1000) ~ Burn_Severity + (1 | Parent_ID)
## model: log10(Leachate_P_mg_kgchar_unfilt/1000) ~ Burn_Severity * Land_Coverage_Category + (1 | Parent_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## model.LC 6 40.336 52.595 -14.1681 28.3363
## model 9 24.615 43.003 -3.3077 6.6153 21.721 3 7.456e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(model,model.int)## refitting model(s) with ML (instead of REML)## Data: data
## Models:
## model.int: log10(Leachate_P_mg_kgchar_unfilt/1000) ~ Burn_Severity + Land_Coverage_Category + (1 | Parent_ID)
## model: log10(Leachate_P_mg_kgchar_unfilt/1000) ~ Burn_Severity * Land_Coverage_Category + (1 | Parent_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## model.int 7 26.705 41.006 -6.3524 12.7048
## model 9 24.615 43.003 -3.3077 6.6153 6.0895 2 0.04761 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#check VIF
vif(model)## GVIF Df GVIF^(1/(2*Df))
## Burn_Severity 2.578947 3 1.171045
## Land_Coverage_Category 5.526316 1 2.350812
## Burn_Severity:Land_Coverage_Category 9.105263 2 1.737093#model assumptions
hist(log10(data$Leachate_P_mg_kgchar_unfilt/1000))
#linearity
ggqqplot(residuals(model))
#check normality of residuals
shapiro.test(log10(data$Leachate_P_mg_kgchar_unfilt)) #p > 0.05 so can assume normality ##
## Shapiro-Wilk normality test
##
## data: log10(data$Leachate_P_mg_kgchar_unfilt)
## W = 0.96458, p-value = 0.09345#check homogeneity of variance assumption (homoscedastcity)
plot(fitted(model),residuals(model))
performance::check_heteroscedasticity(model) #p > 0.05 so can assume homoscedastic## OK: Error variance appears to be homoscedastic (p = 0.913).#post hoc tests for interaction term
emmeans::lsmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Burn_Severity="Unburned"),adjust="tukey") ## $lsmeans
## Burn_Severity Land_Coverage_Category lsmean SE df lower.CL upper.CL
## Unburned Douglas-fir forest -0.885 0.231 12 -1.39 -0.381
## Unburned Sagebrush shrubland -0.501 0.327 12 -1.21 0.211
##
## Degrees-of-freedom method: kenward-roger
## Results are given on the log10 (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df
## (Unburned Douglas-fir forest) - Unburned Sagebrush shrubland -0.384 0.4 12
## t.ratio p.value
## -0.958 0.3571
##
## Degrees-of-freedom method: kenward-roger
## Results are given on the log10 (not the response) scale.emmeans::lsmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Burn_Severity="Low"),adjust="tukey")## $lsmeans
## Burn_Severity Land_Coverage_Category lsmean SE df lower.CL upper.CL
## Low Douglas-fir forest -0.843 0.146 12 -1.162 -0.524
## Low Sagebrush shrubland -0.223 0.231 12 -0.726 0.281
##
## Degrees-of-freedom method: kenward-roger
## Results are given on the log10 (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio
## (Low Douglas-fir forest) - Low Sagebrush shrubland -0.62 0.274 12 -2.268
## p.value
## 0.0426
##
## Degrees-of-freedom method: kenward-roger
## Results are given on the log10 (not the response) scale.emmeans::lsmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Burn_Severity="Moderate"),adjust="tukey")## $lsmeans
## Burn_Severity Land_Coverage_Category lsmean SE df lower.CL upper.CL
## Moderate Douglas-fir forest -0.488 0.189 12 -0.899 -0.0765
## Moderate Sagebrush shrubland 0.833 0.231 12 0.329 1.3368
##
## Degrees-of-freedom method: kenward-roger
## Results are given on the log10 (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df
## (Moderate Douglas-fir forest) - Moderate Sagebrush shrubland -1.32 0.298 12
## t.ratio p.value
## -4.425 0.0008
##
## Degrees-of-freedom method: kenward-roger
## Results are given on the log10 (not the response) scale.emmeans::lsmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Land_Coverage_Category="Douglas-fir forest"),adjust="tukey")## $lsmeans
## Burn_Severity Land_Coverage_Category lsmean SE df lower.CL upper.CL
## High Douglas-fir forest 0.299 0.163 12 -0.0575 0.6549
## Low Douglas-fir forest -0.843 0.146 12 -1.1616 -0.5244
## Moderate Douglas-fir forest -0.488 0.189 12 -0.8991 -0.0765
## Unburned Douglas-fir forest -0.885 0.231 12 -1.3885 -0.3810
##
## Degrees-of-freedom method: kenward-roger
## Results are given on the log10 (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE
## (High Douglas-fir forest) - (Low Douglas-fir forest) 1.1417 0.219
## (High Douglas-fir forest) - (Moderate Douglas-fir forest) 0.7864 0.250
## (High Douglas-fir forest) - (Unburned Douglas-fir forest) 1.1834 0.283
## (Low Douglas-fir forest) - (Moderate Douglas-fir forest) -0.3552 0.239
## (Low Douglas-fir forest) - (Unburned Douglas-fir forest) 0.0418 0.274
## (Moderate Douglas-fir forest) - (Unburned Douglas-fir forest) 0.3970 0.298
## df t.ratio p.value
## 12 5.205 0.0011
## 12 3.149 0.0366
## 12 4.179 0.0061
## 12 -1.488 0.4739
## 12 0.153 0.9987
## 12 1.330 0.5627
##
## Degrees-of-freedom method: kenward-roger
## Results are given on the log10 (not the response) scale.
## P value adjustment: tukey method for comparing a family of 4 estimatesemmeans::lsmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Land_Coverage_Category="Sagebrush shrubland"),adjust="tukey")## $lsmeans
## Burn_Severity Land_Coverage_Category lsmean SE df lower.CL upper.CL
## High Sagebrush shrubland nonEst NA NA NA NA
## Low Sagebrush shrubland -0.223 0.231 12 -0.726 0.281
## Moderate Sagebrush shrubland 0.833 0.231 12 0.329 1.337
## Unburned Sagebrush shrubland -0.501 0.327 12 -1.214 0.211
##
## Degrees-of-freedom method: kenward-roger
## Results are given on the log10 (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df
## High Sagebrush shrubland - Low Sagebrush shrubland nonEst NA NA
## High Sagebrush shrubland - Moderate Sagebrush shrubland nonEst NA NA
## High Sagebrush shrubland - Unburned Sagebrush shrubland nonEst NA NA
## Low Sagebrush shrubland - Moderate Sagebrush shrubland -1.056 0.327 12
## Low Sagebrush shrubland - Unburned Sagebrush shrubland 0.279 0.400 12
## Moderate Sagebrush shrubland - Unburned Sagebrush shrubland 1.334 0.400 12
## t.ratio p.value
## NA NA
## NA NA
## NA NA
## -3.228 0.0318
## 0.696 0.8967
## 3.332 0.0266
##
## Degrees-of-freedom method: kenward-roger
## Results are given on the log10 (not the response) scale.
## P value adjustment: tukey method for comparing a family of 4 estimates#ICP filtered (<0.7 nominal phase) total P stats; likely use this one for publication
model<-lmer(log10(Leachate_P_mg_kgchar_filt0.7/1000)~Burn_Severity*Land_Coverage_Category+(1|Parent_ID), data=data)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficientmodel.BS<-lmer(log10(Leachate_P_mg_kgchar_filt0.7/1000)~Land_Coverage_Category+(1|Parent_ID), data=data)
model.LC<-lmer(log10(Leachate_P_mg_kgchar_filt0.7/1000)~Burn_Severity+(1|Parent_ID), data=data)
model.int<-lmer(log10(Leachate_P_mg_kgchar_filt0.7/1000)~Burn_Severity+Land_Coverage_Category+(1|Parent_ID), data=data)
anova(model,model.BS) #not sig## refitting model(s) with ML (instead of REML)## Data: data
## Models:
## model.BS: log10(Leachate_P_mg_kgchar_filt0.7/1000) ~ Land_Coverage_Category + (1 | Parent_ID)
## model: log10(Leachate_P_mg_kgchar_filt0.7/1000) ~ Burn_Severity * Land_Coverage_Category + (1 | Parent_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## model.BS 4 -37.826 -29.654 22.913 -45.826
## model 9 -33.886 -15.498 25.943 -51.886 6.0598 5 0.3004anova(model,model.LC) #sig## refitting model(s) with ML (instead of REML)## Data: data
## Models:
## model.LC: log10(Leachate_P_mg_kgchar_filt0.7/1000) ~ Burn_Severity + (1 | Parent_ID)
## model: log10(Leachate_P_mg_kgchar_filt0.7/1000) ~ Burn_Severity * Land_Coverage_Category + (1 | Parent_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## model.LC 6 -26.659 -14.401 19.330 -38.659
## model 9 -33.886 -15.498 25.943 -51.886 13.226 3 0.004172 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(model,model.int) #not sig but since BS also isn't sig I probably don't care that much about the LC main effect (so not being able to run the post hoc isn't an issue)## refitting model(s) with ML (instead of REML)## Data: data
## Models:
## model.int: log10(Leachate_P_mg_kgchar_filt0.7/1000) ~ Burn_Severity + Land_Coverage_Category + (1 | Parent_ID)
## model: log10(Leachate_P_mg_kgchar_filt0.7/1000) ~ Burn_Severity * Land_Coverage_Category + (1 | Parent_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## model.int 7 -32.512 -18.211 23.256 -46.512
## model 9 -33.886 -15.498 25.943 -51.886 5.3737 2 0.06809 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#check VIF
vif(model)## GVIF Df GVIF^(1/(2*Df))
## Burn_Severity 2.578947 3 1.171045
## Land_Coverage_Category 5.526316 1 2.350812
## Burn_Severity:Land_Coverage_Category 9.105263 2 1.737093#model assumptions
hist(log10(data$Leachate_P_mg_kgchar_filt0.7/1000))
ggqqplot(residuals(model))
#check normality of residuals
shapiro.test(log10(data$Leachate_P_mg_kgchar_unfilt)) #p > 0.05 so can assume normality ##
## Shapiro-Wilk normality test
##
## data: log10(data$Leachate_P_mg_kgchar_unfilt)
## W = 0.96458, p-value = 0.09345#check homogeneity of variance assumption (homoscedastcity)
plot(fitted(model),residuals(model))
performance::check_heteroscedasticity(model) #need to look into this (fix)## Warning in sqrt(insight::get_deviance(x)/(insight::n_obs(x) -
## sum(!is.na(estimates)))): NaNs produced## Warning: Heteroscedasticity (non-constant error variance) detected (p < .001).#ICP particulate (>0.7 nominal phase) total P stats; likely use this one for publication
model<-lmer(log10(Leachate_P_mg_kgchar_part/1000)~Burn_Severity*Land_Coverage_Category+(1|Parent_ID), data=data)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficientmodel.BS<-lmer(log10(Leachate_P_mg_kgchar_part/1000)~Land_Coverage_Category+(1|Parent_ID), data=data)
model.LC<-lmer(log10(Leachate_P_mg_kgchar_part/1000)~Burn_Severity+(1|Parent_ID), data=data)
model.int<-lmer(log10(Leachate_P_mg_kgchar_part/1000)~Burn_Severity+Land_Coverage_Category+(1|Parent_ID), data=data)
anova(model,model.BS) #sig## refitting model(s) with ML (instead of REML)## Data: data
## Models:
## model.BS: log10(Leachate_P_mg_kgchar_part/1000) ~ Land_Coverage_Category + (1 | Parent_ID)
## model: log10(Leachate_P_mg_kgchar_part/1000) ~ Burn_Severity * Land_Coverage_Category + (1 | Parent_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## model.BS 4 116.60 124.77 -54.298 108.60
## model 9 90.37 108.76 -36.185 72.37 36.226 5 8.559e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(model,model.LC) #sig## refitting model(s) with ML (instead of REML)## Data: data
## Models:
## model.LC: log10(Leachate_P_mg_kgchar_part/1000) ~ Burn_Severity + (1 | Parent_ID)
## model: log10(Leachate_P_mg_kgchar_part/1000) ~ Burn_Severity * Land_Coverage_Category + (1 | Parent_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## model.LC 6 101.81 114.07 -44.905 89.811
## model 9 90.37 108.76 -36.185 72.370 17.441 3 0.0005736 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(model,model.int) #sig## refitting model(s) with ML (instead of REML)## Data: data
## Models:
## model.int: log10(Leachate_P_mg_kgchar_part/1000) ~ Burn_Severity + Land_Coverage_Category + (1 | Parent_ID)
## model: log10(Leachate_P_mg_kgchar_part/1000) ~ Burn_Severity * Land_Coverage_Category + (1 | Parent_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## model.int 7 98.75 113.05 -42.375 84.75
## model 9 90.37 108.76 -36.185 72.37 12.38 2 0.00205 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#check VIF
vif(model)## GVIF Df GVIF^(1/(2*Df))
## Burn_Severity 2.578947 3 1.171045
## Land_Coverage_Category 5.526316 1 2.350812
## Burn_Severity:Land_Coverage_Category 9.105263 2 1.737093#model assumptions
hist(log10(data$Leachate_P_mg_kgchar_filt0.7/1000))
ggqqplot(residuals(model))
#check normality of residuals
shapiro.test(log10(data$Leachate_P_mg_kgchar_unfilt)) #p > 0.05 so can assume normality ##
## Shapiro-Wilk normality test
##
## data: log10(data$Leachate_P_mg_kgchar_unfilt)
## W = 0.96458, p-value = 0.09345#check homogeneity of variance assumption (homoscedastcity)
plot(fitted(model),residuals(model))
performance::check_heteroscedasticity(model) #p > 0.05 so can assume homoscedastic## OK: Error variance appears to be homoscedastic (p = 0.745).#interaction post hoc:
emmeans::lsmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Burn_Severity="Unburned"),adjust="tukey") ## $lsmeans
## Burn_Severity Land_Coverage_Category lsmean SE df lower.CL upper.CL
## Unburned Douglas-fir forest -1.46 0.310 12 -2.14 -0.788
## Unburned Sagebrush shrubland -2.16 0.438 12 -3.11 -1.203
##
## Degrees-of-freedom method: kenward-roger
## Results are given on the log10 (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df
## (Unburned Douglas-fir forest) - Unburned Sagebrush shrubland 0.696 0.537 12
## t.ratio p.value
## 1.296 0.2195
##
## Degrees-of-freedom method: kenward-roger
## Results are given on the log10 (not the response) scale.emmeans::lsmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Burn_Severity="Low"),adjust="tukey") #between ## $lsmeans
## Burn_Severity Land_Coverage_Category lsmean SE df lower.CL upper.CL
## Low Douglas-fir forest -1.322 0.196 12 -1.75 -0.895
## Low Sagebrush shrubland -0.791 0.310 12 -1.47 -0.115
##
## Degrees-of-freedom method: kenward-roger
## Results are given on the log10 (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio
## (Low Douglas-fir forest) - Low Sagebrush shrubland -0.531 0.367 12 -1.449
## p.value
## 0.1731
##
## Degrees-of-freedom method: kenward-roger
## Results are given on the log10 (not the response) scale.emmeans::lsmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Burn_Severity="Moderate"),adjust="tukey") ## $lsmeans
## Burn_Severity Land_Coverage_Category lsmean SE df lower.CL upper.CL
## Moderate Douglas-fir forest -0.687 0.253 12 -1.238 -0.136
## Moderate Sagebrush shrubland 0.825 0.310 12 0.149 1.500
##
## Degrees-of-freedom method: kenward-roger
## Results are given on the log10 (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df
## (Moderate Douglas-fir forest) - Moderate Sagebrush shrubland -1.51 0.4 12
## t.ratio p.value
## -3.778 0.0026
##
## Degrees-of-freedom method: kenward-roger
## Results are given on the log10 (not the response) scale.emmeans::lsmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Land_Coverage_Category="Douglas-fir forest"),adjust="tukey") ## $lsmeans
## Burn_Severity Land_Coverage_Category lsmean SE df lower.CL upper.CL
## High Douglas-fir forest 0.269 0.219 12 -0.208 0.747
## Low Douglas-fir forest -1.322 0.196 12 -1.749 -0.895
## Moderate Douglas-fir forest -0.687 0.253 12 -1.238 -0.136
## Unburned Douglas-fir forest -1.463 0.310 12 -2.138 -0.788
##
## Degrees-of-freedom method: kenward-roger
## Results are given on the log10 (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE
## (High Douglas-fir forest) - (Low Douglas-fir forest) 1.591 0.294
## (High Douglas-fir forest) - (Moderate Douglas-fir forest) 0.956 0.335
## (High Douglas-fir forest) - (Unburned Douglas-fir forest) 1.732 0.380
## (Low Douglas-fir forest) - (Moderate Douglas-fir forest) -0.635 0.320
## (Low Douglas-fir forest) - (Unburned Douglas-fir forest) 0.141 0.367
## (Moderate Douglas-fir forest) - (Unburned Douglas-fir forest) 0.776 0.400
## df t.ratio p.value
## 12 5.410 0.0008
## 12 2.856 0.0608
## 12 4.563 0.0031
## 12 -1.983 0.2472
## 12 0.385 0.9797
## 12 1.939 0.2633
##
## Degrees-of-freedom method: kenward-roger
## Results are given on the log10 (not the response) scale.
## P value adjustment: tukey method for comparing a family of 4 estimatesemmeans::lsmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Land_Coverage_Category="Sagebrush shrubland"),adjust="tukey") ## $lsmeans
## Burn_Severity Land_Coverage_Category lsmean SE df lower.CL upper.CL
## High Sagebrush shrubland nonEst NA NA NA NA
## Low Sagebrush shrubland -0.791 0.310 12 -1.466 -0.115
## Moderate Sagebrush shrubland 0.825 0.310 12 0.149 1.500
## Unburned Sagebrush shrubland -2.159 0.438 12 -3.114 -1.203
##
## Degrees-of-freedom method: kenward-roger
## Results are given on the log10 (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df
## High Sagebrush shrubland - Low Sagebrush shrubland nonEst NA NA
## High Sagebrush shrubland - Moderate Sagebrush shrubland nonEst NA NA
## High Sagebrush shrubland - Unburned Sagebrush shrubland nonEst NA NA
## Low Sagebrush shrubland - Moderate Sagebrush shrubland -1.62 0.438 12
## Low Sagebrush shrubland - Unburned Sagebrush shrubland 1.37 0.537 12
## Moderate Sagebrush shrubland - Unburned Sagebrush shrubland 2.98 0.537 12
## t.ratio p.value
## NA NA
## NA NA
## NA NA
## -3.685 0.0143
## 2.548 0.1020
## 5.557 0.0006
##
## Degrees-of-freedom method: kenward-roger
## Results are given on the log10 (not the response) scale.
## P value adjustment: tukey method for comparing a family of 4 estimates####Solid ICP
solid.P.BS.fig<-ggplot(data_solids, aes(x=Burn_Severity, y=Solid_P_mg_kg, fill=Land_Coverage_Category))+geom_boxplot(show.legend=FALSE)+theme_classic()+ylab(expression(paste(Solid~Total~P~(mg~P~kg^-1), sep="")))+ xlab("Burn Severity")+scale_fill_manual(values=c("#346648", "#dbb40d"))+facet_grid(.~Land_Coverage_Category)+theme(text = element_text(size = 15))
cowplot::save_plot("../figures/solid.P.BS.fig.pdf", solid.P.BS.fig, ncol = 1, nrow = 1, base_aspect_ratio= 2:1, dpi=72)
#ICP solid total P stats; likely use this one for publication
model<-aov(log10(Solid_P_mg_kg)~Burn_Severity*Land_Coverage_Category, data = data_solids)
summary(model)## Df Sum Sq Mean Sq F value Pr(>F)
## Burn_Severity 3 1.0646 0.3549 16.124 0.000165 ***
## Land_Coverage_Category 1 0.7795 0.7795 35.419 6.7e-05 ***
## Burn_Severity:Land_Coverage_Category 2 0.2741 0.1370 6.227 0.013965 *
## Residuals 12 0.2641 0.0220
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#model assumptions
hist(log10(data_solids$Solid_P_mg_kg))
ggqqplot(residuals(model))
#check normality of residuals
shapiro.test(log10(data_solids$Solid_P_mg_kg)) #p > 0.05 so can assume normality ##
## Shapiro-Wilk normality test
##
## data: log10(data_solids$Solid_P_mg_kg)
## W = 0.94907, p-value = 0.3811#check homogeneity of variance assumption (homoscedastcity)
plot(fitted(model),residuals(model))
performance::check_heteroscedasticity(model) #p > 0.05 so can assume homoscedastic## OK: Error variance appears to be homoscedastic (p = 0.558).#interaction post hoc:
emmeans::emmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Burn_Severity="Unburned"),adjust="tukey") ## $emmeans
## Burn_Severity Land_Coverage_Category emmean SE df lower.CL upper.CL
## Unburned Douglas-fir forest 3.09 0.105 12 2.86 3.32
## Unburned Sagebrush shrubland 3.10 0.148 12 2.78 3.43
##
## Results are given on the log10 (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df
## (Unburned Douglas-fir forest) - Unburned Sagebrush shrubland -0.00881 0.182 12
## t.ratio p.value
## -0.048 0.9621
##
## Results are given on the log10 (not the response) scale.emmeans::emmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Burn_Severity="Low"),adjust="tukey") #between ## $emmeans
## Burn_Severity Land_Coverage_Category emmean SE df lower.CL upper.CL
## Low Douglas-fir forest 3.24 0.0663 12 3.09 3.38
## Low Sagebrush shrubland 3.68 0.1049 12 3.45 3.91
##
## Results are given on the log10 (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio
## (Low Douglas-fir forest) - Low Sagebrush shrubland -0.444 0.124 12 -3.579
## p.value
## 0.0038
##
## Results are given on the log10 (not the response) scale.emmeans::emmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Burn_Severity="Moderate"),adjust="tukey") ## $emmeans
## Burn_Severity Land_Coverage_Category emmean SE df lower.CL upper.CL
## Moderate Douglas-fir forest 3.35 0.0857 12 3.16 3.54
## Moderate Sagebrush shrubland 4.15 0.1049 12 3.92 4.38
##
## Results are given on the log10 (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df
## (Moderate Douglas-fir forest) - Moderate Sagebrush shrubland -0.802 0.135 12
## t.ratio p.value
## -5.922 0.0001
##
## Results are given on the log10 (not the response) scale.emmeans::emmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Land_Coverage_Category="Douglas-fir forest"),adjust="tukey") ## $emmeans
## Burn_Severity Land_Coverage_Category emmean SE df lower.CL upper.CL
## High Douglas-fir forest 3.78 0.0742 12 3.61 3.94
## Low Douglas-fir forest 3.24 0.0663 12 3.09 3.38
## Moderate Douglas-fir forest 3.35 0.0857 12 3.16 3.54
## Unburned Douglas-fir forest 3.09 0.1049 12 2.86 3.32
##
## Results are given on the log10 (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE
## (High Douglas-fir forest) - (Low Douglas-fir forest) 0.537 0.0995
## (High Douglas-fir forest) - (Moderate Douglas-fir forest) 0.424 0.1133
## (High Douglas-fir forest) - (Unburned Douglas-fir forest) 0.682 0.1285
## (Low Douglas-fir forest) - (Moderate Douglas-fir forest) -0.113 0.1083
## (Low Douglas-fir forest) - (Unburned Douglas-fir forest) 0.145 0.1241
## (Moderate Douglas-fir forest) - (Unburned Douglas-fir forest) 0.258 0.1354
## df t.ratio p.value
## 12 5.394 0.0008
## 12 3.739 0.0130
## 12 5.305 0.0009
## 12 -1.044 0.7281
## 12 1.167 0.6578
## 12 1.904 0.2768
##
## Results are given on the log10 (not the response) scale.
## P value adjustment: tukey method for comparing a family of 4 estimatesemmeans::emmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Land_Coverage_Category="Sagebrush shrubland"),adjust="tukey") ## $emmeans
## Burn_Severity Land_Coverage_Category emmean SE df lower.CL upper.CL
## High Sagebrush shrubland nonEst NA NA NA NA
## Low Sagebrush shrubland 3.68 0.105 12 3.45 3.91
## Moderate Sagebrush shrubland 4.15 0.105 12 3.92 4.38
## Unburned Sagebrush shrubland 3.10 0.148 12 2.78 3.43
##
## Results are given on the log10 (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df
## High Sagebrush shrubland - Low Sagebrush shrubland nonEst NA NA
## High Sagebrush shrubland - Moderate Sagebrush shrubland nonEst NA NA
## High Sagebrush shrubland - Unburned Sagebrush shrubland nonEst NA NA
## Low Sagebrush shrubland - Moderate Sagebrush shrubland -0.471 0.148 12
## Low Sagebrush shrubland - Unburned Sagebrush shrubland 0.580 0.182 12
## Moderate Sagebrush shrubland - Unburned Sagebrush shrubland 1.051 0.182 12
## t.ratio p.value
## NA NA
## NA NA
## NA NA
## -3.174 0.0350
## 3.193 0.0339
## 5.784 0.0004
##
## Results are given on the log10 (not the response) scale.
## P value adjustment: tukey method for comparing a family of 4 estimates####Correlations
#across both feedstocks
#subset only solid elemental concentrations (in mg/kg)
data_solids_corr <- data_solids %>%
select(matches("Solid"))
# Rename columns by extracting everything between "Solid_" and "mg_Kg"
data_solids_corr <- data_solids_corr %>%
rename_with(
~sub("Solid_(.*?)_mg_kg", "\\1", .),
starts_with("Solid_")) %>%
select(P, C, N, Ca, Mg, Na, Fe, Al, S, K)
#comput r=the correlation matrix, n=the matrix of the number of observations used in analyzing each pari fo variables, and P=the p-values corresponding to the significance levels of correlations
cor_5 <- rcorr(as.matrix(data_solids_corr), type=c("spearman"))
#define M as the correlation coefficient
M <- cor_5$r
#define p_mat as the p-value
p_mat <- cor_5$P
#plot only significant pearson correlation coeffcients when p value >0.05
corrplot(M, type = "upper", p.mat = p_mat, sig.level = 0.05, insig = "blank", tl.col = "black", diag=FALSE, na.label="NA", tl.cex=1, cl.cex = 1, number.cex = 0.8, addCoef.col = "white", col = colorRampPalette(c("#346648", "#dbb40d"))(100))
#Douglas-fir forest
#subset only solid elemental concentrations (in mg/kg)
data_solids_corr <- data_solids %>%
filter(str_detect(Land_Coverage_Category, "Douglas")) %>%
select(matches("Solid"))
# Rename columns by extracting everything between "Solid_" and "mg_Kg"
data_solids_corr <- data_solids_corr %>%
rename_with(
~sub("Solid_(.*?)_mg_kg", "\\1", .),
starts_with("Solid_")) %>%
select(P, C, N, Ca, Mg, Na, Fe, Al, S, K)
#comput r=the correlation matrix, n=the matrix of the number of observations used in analyzing each pari fo variables, and P=the p-values corresponding to the significance levels of correlations
cor_5 <- rcorr(as.matrix(data_solids_corr), type=c("spearman"))
#define M as the correlation coefficient
M <- cor_5$r
#define p_mat as the p-value
p_mat <- cor_5$P
#plot only significant pearson correlation coeffcients when p value >0.05
corrplot(M, type = "upper", p.mat = p_mat, sig.level = 0.05, insig = "blank", tl.col = "black", diag=FALSE, na.label="NA", tl.cex=1, cl.cex = 1, number.cex = 0.8, addCoef.col = "white", col = colorRampPalette(c("#346648", "#dbb40d"))(100))
#Big sagebrush
#subset only solid elemental concentrations (in mg/kg)
data_solids_corr <- data_solids %>%
filter(str_detect(Land_Coverage_Category, "Sagebrush")) %>%
select(matches("Solid"))
# Rename columns by extracting everything between "Solid_" and "mg_Kg"
data_solids_corr <- data_solids_corr %>%
rename_with(
~sub("Solid_(.*?)_mg_kg", "\\1", .),
starts_with("Solid_")) %>%
select(P, C, N, Ca, Mg, Na, Fe, Al, S, K)
#comput r=the correlation matrix, n=the matrix of the number of observations used in analyzing each pari fo variables, and P=the p-values corresponding to the significance levels of correlations
cor_5 <- rcorr(as.matrix(data_solids_corr), type=c("spearman"))
#define M as the correlation coefficient
M <- cor_5$r
#define p_mat as the p-value
p_mat <- cor_5$P
#plot only significant pearson correlation coeffcients when p value >0.05
corrplot(M, type = "upper", p.mat = p_mat, sig.level = 0.05, insig = "blank", tl.col = "black", diag=FALSE, na.label="NA", tl.cex=1, cl.cex = 1, number.cex = 0.8, addCoef.col = "white", col = colorRampPalette(c("#346648", "#dbb40d"))(100))
#P and C
ggplot(data_solids, aes(x=Solid_C_mg_kg, y=Solid_P_mg_kg, color = Land_Coverage_Category)) + geom_point() + geom_smooth(method="lm", se = FALSE) + theme_classic() + scale_color_manual(values = c("#346648", "#dbb40d"))## `geom_smooth()` using formula = 'y ~ x'
#all feedstocks
model<-lm(log(Solid_C_mg_kg)~log(Solid_P_mg_kg), data = data_solids)
summary(model)##
## Call:
## lm(formula = log(Solid_C_mg_kg) ~ log(Solid_P_mg_kg), data = data_solids)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.1667 -0.3490 0.1323 0.3302 0.5910
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.5156 1.0777 12.54 5.11e-10 ***
## log(Solid_P_mg_kg) -0.1815 0.1334 -1.36 0.191
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4742 on 17 degrees of freedom
## Multiple R-squared: 0.09818, Adjusted R-squared: 0.04513
## F-statistic: 1.851 on 1 and 17 DF, p-value: 0.1915#Douglas fir only
model<-lm(log(Solid_C_mg_kg)~log(Solid_P_mg_kg), data = subset(data_solids, Land_Coverage_Category=="Douglas-fir forest"))
summary(model)##
## Call:
## lm(formula = log(Solid_C_mg_kg) ~ log(Solid_P_mg_kg), data = subset(data_solids,
## Land_Coverage_Category == "Douglas-fir forest"))
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.12452 -0.24370 0.00926 0.39848 0.64306
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.2279 1.6388 8.682 1.61e-06 ***
## log(Solid_P_mg_kg) -0.2846 0.2089 -1.363 0.198
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5105 on 12 degrees of freedom
## Multiple R-squared: 0.134, Adjusted R-squared: 0.06181
## F-statistic: 1.856 on 1 and 12 DF, p-value: 0.1981#sagebrush
model<-lm(log(Solid_C_mg_kg)~log(Solid_P_mg_kg), data = subset(data_solids, Land_Coverage_Category=="Sagebrush shrubland"))
summary(model)##
## Call:
## lm(formula = log(Solid_C_mg_kg) ~ log(Solid_P_mg_kg), data = subset(data_solids,
## Land_Coverage_Category == "Sagebrush shrubland"))
##
## Residuals:
## 1 2 3 4 5
## -0.18263 0.07222 -0.17107 0.19096 0.09052
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.76017 0.82811 17.824 0.000385 ***
## log(Solid_P_mg_kg) -0.29516 0.09526 -3.099 0.053357 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1936 on 3 degrees of freedom
## Multiple R-squared: 0.7619, Adjusted R-squared: 0.6826
## F-statistic: 9.601 on 1 and 3 DF, p-value: 0.05336#P and N
ggplot(data_solids, aes(x=Solid_N_mg_kg, y=Solid_P_mg_kg, color = Land_Coverage_Category)) + geom_point() + geom_smooth(method="lm", se = FALSE) + theme_classic() + scale_color_manual(values = c("#346648", "#dbb40d"))## `geom_smooth()` using formula = 'y ~ x'
#all feedstocks
model<-lm(log(Solid_N_mg_kg)~log(Solid_P_mg_kg), data = data_solids)
summary(model)##
## Call:
## lm(formula = log(Solid_N_mg_kg) ~ log(Solid_P_mg_kg), data = data_solids)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.84610 -0.30184 0.02675 0.34999 0.78060
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.08078 1.06949 8.491 1.61e-07 ***
## log(Solid_P_mg_kg) -0.05138 0.13241 -0.388 0.703
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4706 on 17 degrees of freedom
## Multiple R-squared: 0.008781, Adjusted R-squared: -0.04953
## F-statistic: 0.1506 on 1 and 17 DF, p-value: 0.7028#Douglas fir only
model<-lm(log(Solid_N_mg_kg)~log(Solid_P_mg_kg), data = subset(data_solids, Land_Coverage_Category=="Douglas-fir forest"))
summary(model)##
## Call:
## lm(formula = log(Solid_N_mg_kg) ~ log(Solid_P_mg_kg), data = subset(data_solids,
## Land_Coverage_Category == "Douglas-fir forest"))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.78941 -0.30421 -0.04207 0.37311 0.84806
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.4408 1.6029 5.890 7.37e-05 ***
## log(Solid_P_mg_kg) -0.1090 0.2043 -0.534 0.603
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4993 on 12 degrees of freedom
## Multiple R-squared: 0.02318, Adjusted R-squared: -0.05822
## F-statistic: 0.2848 on 1 and 12 DF, p-value: 0.6033#sagebrush
model<-lm(log(Solid_N_mg_kg)~log(Solid_P_mg_kg), data = subset(data_solids, Land_Coverage_Category=="Sagebrush shrubland"))
summary(model)##
## Call:
## lm(formula = log(Solid_N_mg_kg) ~ log(Solid_P_mg_kg), data = subset(data_solids,
## Land_Coverage_Category == "Sagebrush shrubland"))
##
## Residuals:
## 1 2 3 4 5
## -0.02507 0.31641 0.11150 -0.21428 -0.18856
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.8132 1.0893 9.927 0.00217 **
## log(Solid_P_mg_kg) -0.2224 0.1253 -1.775 0.17400
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2547 on 3 degrees of freedom
## Multiple R-squared: 0.5122, Adjusted R-squared: 0.3496
## F-statistic: 3.15 on 1 and 3 DF, p-value: 0.174####PCA & RDAs
#PCA of elemental concentration of the solid chars with burn severity and feedstock
#format dataframe using same code as correlation and then also remove NAs
data_solids_corr <- data_solids %>%
rename_with(
~sub("Solid_(.*?)_mg_kg", "\\1", .),
starts_with("Solid_")) %>%
select(Burn_Severity, Land_Coverage_Category, P, C, N, Ca, Mg, Na, Fe, Al, S, K) %>%
na.omit(data_solids)
pca <- prcomp(data_solids_corr[3:12], center=TRUE, scale=TRUE)
##Extract PCA information, including loadings, sample scores, and component percentages
eigval <- pca$sdev^2
eigvec <- data.frame(pca$rotation)
loadings1 <- eigvec$PC1*sqrt(eigval[1])
loadings2 <- eigvec$PC2*sqrt(eigval[2])
loadings <- cbind(loadings1, loadings2)
colnames(loadings)<-c('PC1','PC2')
rownames(loadings)<-colnames(data_solids_corr[3:12])
loadingstest<-as.data.frame(loadings)
pca_out <- as.data.frame(pca$x)
percentage <- round(pca$sdev^2/sum(pca$sdev^2)*100,1)
percentage <- paste(colnames(pca_out), "-", paste(as.character(percentage), "%", sep=""))
xscale <- 10
yscale <- 10
ggplot()+geom_point(data=pca_out,mapping=aes(x=PC1, y=PC2, colour=data_solids_corr$Land_Coverage_Category,shape=data_solids_corr$Burn_Severity),size=3)+xlim(-xscale,xscale)+ylim(-yscale,yscale)+geom_abline(intercept=0, slope=0, linetype="dashed", size=0.8, colour="gray")+geom_vline(aes(xintercept=0), linetype="dashed", size=0.8, colour="gray")+theme_light()+labs(colour="Land_Coverage_Category",shape="Burn_Severity")+theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank(),panel.background=element_blank(),axis.title.x=element_text(size=16),axis.title.y=element_text(size=16),axis.text.y=element_text(size=16),axis.text.x=element_text(size=16),legend.title=element_text(size=12),legend.text=element_text(size=12),legend.position="bottom",legend.box="verticle",legend.margin=margin())+xlab(percentage[1])+ylab(percentage[2])+guides(color=guide_legend(nrow=2, byrow=TRUE))+scale_color_manual(values=c("#346648", "#dbb40d"))+scale_shape_manual(values=c(16,17,3,12))+geom_text_repel(data=loadingstest,mapping=aes(label=row.names(loadingstest), x=PC1*xscale, y=PC2*yscale),size=5, colour="black", alpha=0.5)
#include pH, particulate concentrations from ICP, and filt0.7 concentrations from ICP, NPOC, TDN, MBD
#fix THIS DOESN"T WORK RIGHT NOW
#format dataframe and then also remove NAs
#data_solids_corr <- data_solids %>%
# rename_with(
# ~sub("Leachate_(.*?)_mg_kg", "\\1", .),
# starts_with("Leachate_")) %>%
# select(Burn_Severity, Land_Coverage_Category, P, C, N, Ca, Mg, Na, Fe, Al, S, K) %>%
# na.omit(data_solids)#PCAs of XANES LCF
#FIX THIS DOENS'T RUN RIGHT NOW
#all phases together; open air; Po individual species
#format dataframe using same code as correlation and then also remove NAs
#XANES_corr <- XANES_clean %>%
# select(Parent_ID, Burn_Severity, Land_Coverage_Category, FilterSize, XANES_Pi_Al, XANES_Pi_Ca, XANES_Pi_Fe, #XANES_Pi_K, XANES_Pi_Mg, XANES_Pi_Na, XANES_Po_Al, XANES_Po_Fe, XANES_Po_Ca, XANES_Po_Na) %>%
# na.omit(XANES_clean)
#pca <- prcomp(XANES_corr[5:14], center=TRUE, scale=TRUE)
##Extract PCA information, including loadings, sample scores, and component percentages
#eigval <- pca$sdev^2
#eigvec <- data.frame(pca$rotation)
#loadings1 <- eigvec$PC1*sqrt(eigval[1])
#loadings2 <- eigvec$PC2*sqrt(eigval[2])
#loadings <- cbind(loadings1, loadings2)
#colnames(loadings)<-c('PC1','PC2')
#rownames(loadings)<-colnames(XANES_corr[5:14])
#loadingstest<-as.data.frame(loadings)
#pca_out <- as.data.frame(pca$x)
#percentage <- round(pca$sdev^2/sum(pca$sdev^2)*100,1)
#percentage <- paste(colnames(pca_out), "-", paste(as.character(percentage), "%", sep=""))
#xscale <- 4
#yscale <- 4
#ggplot()+geom_point(data=pca_out,mapping=aes(x=PC1, y=PC2, colour=XANES_corr$Land_Coverage_Category,shape=XANES_corr$Burn_Severity, size=XANES_corr$FilterSize))+xlim(-xscale,xscale)+ylim(-yscale,yscale)+geom_abline(intercept=0, slope=0, linetype="dashed", size=0.8, colour="gray")+geom_vline(aes(xintercept=0), linetype="dashed", size=0.8, colour="gray")+theme_light()+labs(colour="Land_Coverage_Category",shape="Burn_Severity")+theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank(),panel.background=element_blank(),axis.title.x=element_text(size=16),axis.title.y=element_text(size=16),axis.text.y=element_text(size=16),axis.text.x=element_text(size=16),legend.title=element_text(size=12),legend.text=element_text(size=12),legend.position="bottom",legend.box="verticle",legend.margin=margin())+xlab(percentage[1])+ylab(percentage[2])+guides(color=guide_legend(nrow=2, byrow=TRUE))+scale_color_manual(values=c("#346648", "#dbb40d"))+scale_shape_manual(values=c(16,17,3,12))+scale_size_manual(values=c(3,6,9))+geom_text_repel(data=loadingstest,mapping=aes(label=row.names(loadingstest), x=PC1*xscale, y=PC2*yscale),size=5, colour="black", alpha=0.5)+ggtitle("Open Air XANES - all phases")
#all phases together; open air; Po total instead of individual
#format dataframe using same code as correlation and then also remove NAs
#XANES_corr <- XANES_clean %>%
# select(Parent_ID, Burn_Severity, Land_Coverage_Category, FilterSize, XANES_Pi_Al, XANES_Pi_Ca, XANES_Pi_Fe, XANES_Pi_K, XANES_Pi_Mg, XANES_Pi_Na, XANES_Po) %>%
# na.omit(XANES_clean)
#pca <- prcomp(XANES_corr[5:11], center=TRUE, scale=TRUE)
##Extract PCA information, including loadings, sample scores, and component percentages
#eigval <- pca$sdev^2
#eigvec <- data.frame(pca$rotation)
#loadings1 <- eigvec$PC1*sqrt(eigval[1])
#loadings2 <- eigvec$PC2*sqrt(eigval[2])
#loadings <- cbind(loadings1, loadings2)
#colnames(loadings)<-c('PC1','PC2')
#rownames(loadings)<-colnames(XANES_corr[5:11])
#loadingstest<-as.data.frame(loadings)
#pca_out <- as.data.frame(pca$x)
#percentage <- round(pca$sdev^2/sum(pca$sdev^2)*100,1)
#percentage <- paste(colnames(pca_out), "-", paste(as.character(percentage), "%", sep=""))
#xscale <- 4
#yscale <- 4
#ggplot()+geom_point(data=pca_out,mapping=aes(x=PC1, y=PC2, colour=XANES_corr$Land_Coverage_Category,shape=XANES_corr$Burn_Severity, size=XANES_corr$FilterSize))+xlim(-xscale,xscale)+ylim(-yscale,yscale)+geom_abline(intercept=0, slope=0, linetype="dashed", size=0.8, colour="gray")+geom_vline(aes(xintercept=0), linetype="dashed", size=0.8, colour="gray")+theme_light()+labs(colour="Land_Coverage_Category",shape="Burn_Severity")+theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank(),panel.background=element_blank(),axis.title.x=element_text(size=16),axis.title.y=element_text(size=16),axis.text.y=element_text(size=16),axis.text.x=element_text(size=16),legend.title=element_text(size=12),legend.text=element_text(size=12),legend.position="bottom",legend.box="verticle",legend.margin=margin())+xlab(percentage[1])+ylab(percentage[2])+guides(color=guide_legend(nrow=2, byrow=TRUE))+scale_color_manual(values=c("#346648", "#dbb40d"))+scale_shape_manual(values=c(16,17,3,12))+scale_size_manual(values=c(3,6,9))+geom_text_repel(data=loadingstest,mapping=aes(label=row.names(loadingstest), x=PC1*xscale, y=PC2*yscale),size=5, colour="black", alpha=0.5)+ggtitle("Open Air XANES - all phases")
#solid only; open air; Po individual species
#format dataframe using same code as correlation and then also remove NAs; also removed columns that were all zeros for this phase
XANES_corr <- XANES_clean %>%
select(Parent_ID, Burn_Severity, Land_Coverage_Category, FilterSize, XANES_Pi_Al, XANES_Pi_Ca, XANES_Pi_Fe, XANES_Pi_K, XANES_Pi_Mg, XANES_Pi_Na, XANES_Po_Ca, XANES_Po_Na) %>%
na.omit(XANES_clean) %>%
filter(FilterSize %in% c("solid"))
pca <- prcomp(XANES_corr[5:12], center=TRUE, scale=TRUE)
##Extract PCA information, including loadings, sample scores, and component percentages
eigval <- pca$sdev^2
eigvec <- data.frame(pca$rotation)
loadings1 <- eigvec$PC1*sqrt(eigval[1])
loadings2 <- eigvec$PC2*sqrt(eigval[2])
loadings <- cbind(loadings1, loadings2)
colnames(loadings)<-c('PC1','PC2')
rownames(loadings)<-colnames(XANES_corr[5:12])
loadingstest<-as.data.frame(loadings)
pca_out <- as.data.frame(pca$x)
percentage <- round(pca$sdev^2/sum(pca$sdev^2)*100,1)
percentage <- paste(colnames(pca_out), "-", paste(as.character(percentage), "%", sep=""))
xscale <- 4
yscale <- 4
ggplot()+geom_point(data=pca_out,mapping=aes(x=PC1, y=PC2, colour=XANES_corr$Land_Coverage_Category,shape=XANES_corr$Burn_Severity),size=5)+xlim(-xscale,xscale)+ylim(-yscale,yscale)+geom_abline(intercept=0, slope=0, linetype="dashed", size=0.8, colour="gray")+geom_vline(aes(xintercept=0), linetype="dashed", size=0.8, colour="gray")+theme_light()+labs(colour="Land_Coverage_Category",shape="Burn_Severity")+theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank(),panel.background=element_blank(),axis.title.x=element_text(size=16),axis.title.y=element_text(size=16),axis.text.y=element_text(size=16),axis.text.x=element_text(size=16),legend.title=element_text(size=12),legend.text=element_text(size=12),legend.position="bottom",legend.box="verticle",legend.margin=margin())+xlab(percentage[1])+ylab(percentage[2])+guides(color=guide_legend(nrow=2, byrow=TRUE))+scale_color_manual(values=c("#346648", "#dbb40d"))+scale_shape_manual(values=c(16,17,3,12))+geom_text_repel(data=loadingstest,mapping=aes(label=row.names(loadingstest), x=PC1*xscale, y=PC2*yscale),size=5, colour="black", alpha=0.5)+ggtitle("Open Air XANES - Solid")
#solid only; open air; Po grouped
#format dataframe using same code as correlation and then also remove NAs; also removed columns that were all zeros for this phase
XANES_corr <- XANES_clean %>%
select(Parent_ID, Burn_Severity, Land_Coverage_Category, FilterSize, XANES_Pi_Al, XANES_Pi_Ca, XANES_Pi_Fe, XANES_Pi_K, XANES_Pi_Mg, XANES_Pi_Na, XANES_Po) %>%
na.omit(XANES_clean) %>%
filter(FilterSize %in% c("solid"))
pca <- prcomp(XANES_corr[5:11], center=TRUE, scale=TRUE)
##Extract PCA information, including loadings, sample scores, and component percentages
eigval <- pca$sdev^2
eigvec <- data.frame(pca$rotation)
loadings1 <- eigvec$PC1*sqrt(eigval[1])
loadings2 <- eigvec$PC2*sqrt(eigval[2])
loadings <- cbind(loadings1, loadings2)
colnames(loadings)<-c('PC1','PC2')
rownames(loadings)<-colnames(XANES_corr[5:11])
loadingstest<-as.data.frame(loadings)
pca_out <- as.data.frame(pca$x)
percentage <- round(pca$sdev^2/sum(pca$sdev^2)*100,1)
percentage <- paste(colnames(pca_out), "-", paste(as.character(percentage), "%", sep=""))
xscale <- 4
yscale <- 4
ggplot()+geom_point(data=pca_out,mapping=aes(x=PC1, y=PC2, colour=XANES_corr$Land_Coverage_Category,shape=XANES_corr$Burn_Severity),size=5)+xlim(-xscale,xscale)+ylim(-yscale,yscale)+geom_abline(intercept=0, slope=0, linetype="dashed", size=0.8, colour="gray")+geom_vline(aes(xintercept=0), linetype="dashed", size=0.8, colour="gray")+theme_light()+labs(colour="Land_Coverage_Category",shape="Burn_Severity")+theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank(),panel.background=element_blank(),axis.title.x=element_text(size=16),axis.title.y=element_text(size=16),axis.text.y=element_text(size=16),axis.text.x=element_text(size=16),legend.title=element_text(size=12),legend.text=element_text(size=12),legend.position="bottom",legend.box="verticle",legend.margin=margin())+xlab(percentage[1])+ylab(percentage[2])+guides(color=guide_legend(nrow=2, byrow=TRUE))+scale_color_manual(values=c("#346648", "#dbb40d"))+scale_shape_manual(values=c(16,17,3,12))+geom_text_repel(data=loadingstest,mapping=aes(label=row.names(loadingstest), x=PC1*xscale, y=PC2*yscale),size=5, colour="black", alpha=0.5)+ggtitle("Open Air XANES - Solid")
leachate <- data %>%
select(1:3, 5, 10, 11, 49, 62, 70, 78, 50)
leachate_mean <- leachate %>%
dplyr::group_by(Parent_ID) %>%
dplyr::summarize(
Leachate_P_mg_kgchar_unfilt_Mean = mean(Leachate_P_mg_kgchar_unfilt, na.rm = TRUE),
Leachate_P_mg_kgchar_part_Mean = mean(Leachate_P_mg_kgchar_part, na.rm = TRUE),
Leachate_P_mg_kgchar_filt0.7_Mean = mean(Leachate_P_mg_kgchar_filt0.7, na.rm = TRUE),
MBD_P_mg_kgchar_Mean = mean(MBD_P_mg_kgchar, na.rm = TRUE),
pH_Mean = mean(pH)
)
solid <- data_solids %>%
select(1,3, 5, 10, 11, 31:34, 88, 96:101, 111:120)
data_rda<- solid %>%
full_join(leachate_mean) #%>%## Joining with `by = join_by(Parent_ID)` #na.omit(solid) #remove NAs
data_rda <- data_rda %>%
select(2,3,6:17, 27:30)
pca <- prcomp(data_rda[3:18], center=TRUE, scale=TRUE)
##Extract PCA information, including loadings, sample scores, and component percentages
eigval <- pca$sdev^2
eigvec <- data.frame(pca$rotation)
loadings1 <- eigvec$PC1*sqrt(eigval[1])
loadings2 <- eigvec$PC2*sqrt(eigval[2])
loadings <- cbind(loadings1, loadings2)
colnames(loadings)<-c('PC1','PC2')
rownames(loadings)<-colnames(data_rda[3:18])
loadingstest<-as.data.frame(loadings)
pca_out <- as.data.frame(pca$x)
percentage <- round(pca$sdev^2/sum(pca$sdev^2)*100,1)
percentage <- paste(colnames(pca_out), "-", paste(as.character(percentage), "%", sep=""))
xscale <- 6
yscale <- 6
ggplot()+geom_point(data=pca_out,mapping=aes(x=PC1, y=PC2, colour=data_rda$Land_Coverage_Category,shape=data_rda$Burn_Severity),size=5)+xlim(-xscale,xscale)+ylim(-yscale,yscale)+geom_abline(intercept=0, slope=0, linetype="dashed", size=0.8, colour="gray")+geom_vline(aes(xintercept=0), linetype="dashed", size=0.8, colour="gray")+theme_light()+labs(colour="Land_Coverage_Category",shape="Burn_Severity")+theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank(),panel.background=element_blank(),axis.title.x=element_text(size=16),axis.title.y=element_text(size=16),axis.text.y=element_text(size=16),axis.text.x=element_text(size=16),legend.title=element_text(size=12),legend.text=element_text(size=12),legend.position="bottom",legend.box="verticle",legend.margin=margin())+xlab(percentage[1])+ylab(percentage[2])+guides(color=guide_legend(nrow=2, byrow=TRUE))+scale_color_manual(values=c("#346648", "#dbb40d"))+scale_shape_manual(values=c(16,17,3,12))+geom_text_repel(data=loadingstest,mapping=aes(label=row.names(loadingstest), x=PC1*xscale, y=PC2*yscale),size=5, colour="black", alpha=0.5)
#setting up an RDA to look at how burn conditions (severity, duration, max temp) and chemistry of solid char (P, Ca, Mg, K, Fe, Al, Na concentration, XANES categories, and NMR categories) impact P mobilization (total P concentration in unfiltered leachate, aqueous phase (<0.7), and particulate phase (>0.7); soluble reactive P)
#this is example code of how to do an RDA
#reference: http://dmcglinn.github.io/quant_methods/lessons/multivariate_models.html
#other reference: https://popgen.nescent.org/2018-03-27_RDA_GEA.html
#citation in here: https://cran.r-project.org/web/packages/vegan/vignettes/FAQ-vegan.html
#reference about scaling: https://mb3is.megx.net/gustame/constrained-analyses/rda
#make dataframe with all leachate samples
leachate <- data %>%
select(1:3, 5, 10, 11, 49, 62, 70, 78, 50)
solid <- data_solids %>%
select(1, 31:34, 88, 96:101, 111:120)
data_rda<- leachate %>%
full_join(solid) %>%
filter(Solid_Na_g_kg != "") #remove NAs## Joining with `by = join_by(Parent_ID)`#create dataframe for sample variables
samples<-data_rda[ c(7:32)]
#create dataframe for burn condition variables
enviro<-data_rda[ c(4:6)]
#make dataframe with solid samples (leachate avg) to prevent using same number so many times
#alternative RDA: avergae leachate concentrations
leachate_mean <- leachate %>%
dplyr::group_by(Parent_ID) %>%
dplyr::summarize(
Leachate_P_mg_kgchar_unfilt_Mean = mean(Leachate_P_mg_kgchar_unfilt, na.rm = TRUE),
Leachate_P_mg_kgchar_part_Mean = mean(Leachate_P_mg_kgchar_part, na.rm = TRUE),
Leachate_P_mg_kgchar_filt0.7_Mean = mean(Leachate_P_mg_kgchar_filt0.7, na.rm = TRUE),
MBD_P_mg_kgchar_Mean = mean(MBD_P_mg_kgchar, na.rm = TRUE),
pH_Mean = mean(pH)
)
solid <- data_solids %>%
select(1,3, 5, 10, 11, 31:34, 88, 96:101, 111:120)
data_rda<- solid %>%
full_join(leachate_mean) %>%
filter(Solid_Na_g_kg != "") #remove NAs## Joining with `by = join_by(Parent_ID)`#RDA: How do burn conditions (severity, duration, temp) influence chemistry of solids and leachates?
#create dataframe for sample variables
samples<-data_rda[ c(6:31)]
#create dataframe for burn condition variables
enviro<-data_rda[ c(3:5)]
#make sure both dataframes have equal rows
all.equal(rownames(samples), rownames(enviro))## [1] TRUE#RDA (this anaysis expects a linear response of each sample to the environmental variables)
rda_tree=rda(samples~ . , data=enviro)
rda_tree## Call: rda(formula = samples ~ Burn_Severity + Burn_Duration +
## Char_Max_Temp, data = enviro)
##
## Inertia Proportion Rank
## Total 1.975e+07 1.000e+00
## Constrained 7.086e+06 3.588e-01 5
## Unconstrained 1.267e+07 6.412e-01 12
## Inertia is variance
##
## Eigenvalues for constrained axes:
## RDA1 RDA2 RDA3 RDA4 RDA5
## 7077630 5994 2089 345 87
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
## 12637555 26075 1237 466 212 133 77 59
## PC9 PC10 PC11 PC12
## 20 5 1 1RsquareAdj(rda_tree)## $r.squared
## [1] 0.358756
##
## $adj.r.squared
## [1] 0.09157107summary(eigenvals(rda_tree, model = "constrained"))## Importance of components:
## RDA1 RDA2 RDA3 RDA4 RDA5
## Eigenvalue 7.078e+06 5.994e+03 2.089e+03 3.454e+02 8.651e+01
## Proportion Explained 9.988e-01 8.458e-04 2.948e-04 4.875e-05 1.221e-05
## Cumulative Proportion 9.988e-01 9.996e-01 9.999e-01 1.000e+00 1.000e+00#see the factor loadings (correlations of each variable with each pca axis)
scores(rda_tree, choices=1:2, display="species", scaling=2)## RDA1 RDA2
## Ortho -0.258608414 -0.3836313766
## Mono 0.174726143 0.3077998691
## Di 0.058838394 0.1025543877
## Pyro 0.025383379 -0.0272017518
## XANES_Po 0.281408979 0.5249502154
## XANES_Pi_Fe 0.118469141 0.3483344631
## XANES_Pi_Na 0.209494308 -0.2262644254
## XANES_Pi_Ca -0.248376554 -0.4292863568
## XANES_Pi_Al -0.011770369 -0.0002355077
## XANES_Pi_Mg -0.317177075 -0.1469297923
## XANES_Pi_K -0.014982768 -0.0554861412
## Solid_P_g_kg -0.080029673 -0.0212126139
## Solid_Ca_g_kg -0.365725846 -0.3100954065
## Solid_Fe_g_kg -0.050220016 -0.0196842734
## Solid_Al_g_kg -0.047826939 -0.0208178893
## Solid_K_g_kg -0.558503286 -0.0876518422
## Solid_Mg_g_kg -0.073643548 -0.0213093323
## Solid_Na_g_kg -0.009612414 -0.0046855763
## Solid_S_g_kg -0.021646259 0.0001316653
## Solid_C_g_kg 1.091870534 2.1034905286
## Solid_N_g_kg 0.019330517 0.0613465720
## Leachate_P_mg_kgchar_unfilt_Mean -56.778548158 0.1556707251
## Leachate_P_mg_kgchar_part_Mean -57.781948472 -0.0964623423
## Leachate_P_mg_kgchar_filt0.7_Mean 1.003400315 0.2521330673
## MBD_P_mg_kgchar_Mean 0.684190842 0.1321200187
## pH_Mean -0.050300766 -0.0401592675
## attr(,"const")
## [1] 135.3677#plot it
plot(rda_tree, type='n', scaling=2)
orditorp(rda_tree, display='sp', cex=0.5, scaling=1, col='blue')
text(rda_tree, display='cn', col='red')
#check correlation
pairs.panels(enviro, scale=T)
#can visualize RDA importance
screeplot(rda_tree)
#check significance of all PCAs
signif.full <- anova.cca(rda_tree, parallel=getOption("mc.cores")) # default is permutation=999
signif.full #model is significant## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = samples ~ Burn_Severity + Burn_Duration + Char_Max_Temp, data = enviro)
## Df Variance F Pr(>F)
## Model 5 7086145 1.3427 0.259
## Residual 12 12665842#check significance of individual PCAs
signif.axis <- anova.cca(rda_tree, by="axis", parallel=getOption("mc.cores"))
signif.axis ## Permutation test for rda under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = samples ~ Burn_Severity + Burn_Duration + Char_Max_Temp, data = enviro)
## Df Variance F Pr(>F)
## RDA1 1 7077630 6.7056 0.272
## RDA2 1 5994 0.0057 1.000
## RDA3 1 2089 0.0020 1.000
## RDA4 1 345 0.0003 1.000
## RDA5 1 87 0.0001 1.000
## Residual 12 12665842#check VIF
vif.cca(rda_tree)## Burn_SeverityLow Burn_SeverityModerate Burn_SeverityUnburned
## 3.437847 2.249385 8.802287
## Burn_Duration Char_Max_Temp
## 2.638543 8.735625#compare the model to chance with permutation test
anova(rda_tree)## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = samples ~ Burn_Severity + Burn_Duration + Char_Max_Temp, data = enviro)
## Df Variance F Pr(>F)
## Model 5 7086145 1.3427 0.28
## Residual 12 12665842#plot again
plot(rda_tree, scaling=2) #this just doesn't work; not significant, VIF too high
#RDA: How does solid P chemistry (NMR, XANES) influence elemental concentrations in leachate (unfilt, filt, part, molybdate)?; color by burn severity, shape by feedstock; uses leachate means; VIF is too high
data_rda<- solid %>%
full_join(leachate_mean) %>%
filter(Solid_Na_g_kg != "") #remove NAs## Joining with `by = join_by(Parent_ID)`#create dataframe for sample variables
samples<-data_rda[ c(27:30)]
#create dataframe for burn condition variables
enviro<-data_rda[ c(6:16)]
#make sure both dataframes have equal rows
all.equal(rownames(samples), rownames(enviro))## [1] TRUE#RDA (this anaysis expects a linear response of each sample to the environmental variables)
rda_tree=rda(samples~ . , data=enviro)
rda_tree## Call: rda(formula = samples ~ Ortho + Mono + Di + Pyro + XANES_Po +
## XANES_Pi_Fe + XANES_Pi_Na + XANES_Pi_Ca + XANES_Pi_Al + XANES_Pi_Mg +
## XANES_Pi_K, data = enviro)
##
## Inertia Proportion Rank
## Total 1.974e+07 1.000e+00
## Constrained 1.651e+07 8.366e-01 3
## Unconstrained 3.225e+06 1.634e-01 3
## Inertia is variance
##
## Eigenvalues for constrained axes:
## RDA1 RDA2 RDA3
## 16497551 16809 359
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3
## 3213434 12019 38RsquareAdj(rda_tree)## $r.squared
## [1] 0.836603
##
## $adj.r.squared
## [1] 0.5370419summary(eigenvals(rda_tree, model = "constrained"))## Importance of components:
## RDA1 RDA2 RDA3
## Eigenvalue 1.65e+07 1.681e+04 3.591e+02
## Proportion Explained 9.99e-01 1.018e-03 2.174e-05
## Cumulative Proportion 9.99e-01 1.000e+00 1.000e+00#see the factor loadings (correlations of each variable with each pca axis)
scores(rda_tree, choices=1:2, display="species", scaling=2)## RDA1 RDA2
## Leachate_P_mg_kgchar_unfilt_Mean 87.3227700 -1.432921
## Leachate_P_mg_kgchar_part_Mean 87.6604394 1.413790
## Leachate_P_mg_kgchar_filt0.7_Mean -0.3376693 -2.846711
## MBD_P_mg_kgchar_Mean -0.1250006 -1.855479
## attr(,"const")
## [1] 135.3475#plot it
plot(rda_tree, type='n', scaling=2)
orditorp(rda_tree, display='sp', cex=0.5, scaling=1, col='blue')
text(rda_tree, display='cn', col='red')
#check correlation
pairs.panels(enviro, scale=T)
#can visualize RDA importance
screeplot(rda_tree)
#check significance of all PCAs
signif.full <- anova.cca(rda_tree, parallel=getOption("mc.cores")) # default is permutation=999
signif.full #model is significant## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = samples ~ Ortho + Mono + Di + Pyro + XANES_Po + XANES_Pi_Fe + XANES_Pi_Na + XANES_Pi_Ca + XANES_Pi_Al + XANES_Pi_Mg + XANES_Pi_K, data = enviro)
## Df Variance F Pr(>F)
## Model 11 16514719 2.7928 0.133
## Residual 6 3225491#check significance of individual PCAs
signif.axis <- anova.cca(rda_tree, by="axis", parallel=getOption("mc.cores"))
signif.axis #only RDA1 and RDA 2 are sig## Permutation test for rda under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = samples ~ Ortho + Mono + Di + Pyro + XANES_Po + XANES_Pi_Fe + XANES_Pi_Na + XANES_Pi_Ca + XANES_Pi_Al + XANES_Pi_Mg + XANES_Pi_K, data = enviro)
## Df Variance F Pr(>F)
## RDA1 1 16497551 71.6064 0.121
## RDA2 1 16809 0.0730 1.000
## RDA3 1 359 0.0016 1.000
## Residual 14 3225491#check VIF
vif.cca(rda_tree)## Ortho Mono Di Pyro XANES_Po XANES_Pi_Fe
## 430897.99439 244137.82681 27453.58296 31271.68532 2066.14942 1017.06394
## XANES_Pi_Na XANES_Pi_Ca XANES_Pi_Al XANES_Pi_Mg XANES_Pi_K
## 2069.52735 1022.25453 83.70275 1425.22502 481.65760#compare the model to chance with permutation test
anova(rda_tree)## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = samples ~ Ortho + Mono + Di + Pyro + XANES_Po + XANES_Pi_Fe + XANES_Pi_Na + XANES_Pi_Ca + XANES_Pi_Al + XANES_Pi_Mg + XANES_Pi_K, data = enviro)
## Df Variance F Pr(>F)
## Model 11 16514719 2.7928 0.135
## Residual 6 3225491#plot again
plot(rda_tree, scaling=2) #this just doesn't work; not significant, VIF too high
#How does feedstock and burn severity influence solid P chemistry (P conc, NMR, XANES)?
data_rda<- solid %>%
full_join(leachate_mean) ## Joining with `by = join_by(Parent_ID)`#create dataframe for sample variables
samples<-data_rda[ c(6:17)]
#create dataframe for burn condition variables
enviro<-data_rda[ c(2:3)]
#make sure both dataframes have equal rows
all.equal(rownames(samples), rownames(enviro))## [1] TRUE#RDA (this anaysis expects a linear response of each sample to the environmental variables)
rda_tree=rda(samples~ . , data=enviro)
rda_tree## Call: rda(formula = samples ~ Land_Coverage_Category + Burn_Severity,
## data = enviro)
##
## Inertia Proportion Rank
## Total 2544.1397 1.0000
## Constrained 1767.1852 0.6946 4
## Unconstrained 776.9546 0.3054 12
## Inertia is variance
##
## Eigenvalues for constrained axes:
## RDA1 RDA2 RDA3 RDA4
## 1300.3 399.3 37.0 30.6
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10 PC11 PC12
## 420.4 147.8 83.8 71.1 31.6 16.2 2.9 2.3 0.7 0.3 0.0 0.0RsquareAdj(rda_tree)## $r.squared
## [1] 0.6946101
##
## $adj.r.squared
## [1] 0.6073559summary(eigenvals(rda_tree, model = "constrained"))## Importance of components:
## RDA1 RDA2 RDA3 RDA4
## Eigenvalue 1300.2585 399.2663 37.02949 30.63083
## Proportion Explained 0.7358 0.2259 0.02095 0.01733
## Cumulative Proportion 0.7358 0.9617 0.98267 1.00000#see the factor loadings (correlations of each variable with each pca axis)
scores(rda_tree, choices=1:2, display="species", scaling=2)## RDA1 RDA2
## Ortho 4.12832420 1.1490159
## Mono -3.24044129 -0.2501559
## Di -1.20235534 0.2376141
## Pyro 0.31901678 -1.1398120
## XANES_Po -5.78118619 0.9835438
## XANES_Pi_Fe -3.73794387 1.3954752
## XANES_Pi_Na 0.88496193 -4.8047389
## XANES_Pi_Ca 4.49340039 0.7683397
## XANES_Pi_Al 0.08413426 -0.4895688
## XANES_Pi_Mg 3.35883295 1.8850168
## XANES_Pi_K 0.50231586 0.2399660
## Solid_P_g_kg 0.60750473 0.6357070
## attr(,"const")
## [1] 14.62862#plot it
plot(rda_tree, type='n', scaling=2)
orditorp(rda_tree, display='sp', cex=0.5, scaling=1, col='blue')
text(rda_tree, display='cn', col='red')
#check correlation
pairs.panels(enviro, scale=T)
#can visualize RDA importance
screeplot(rda_tree)
#check significance of all PCAs
signif.full <- anova.cca(rda_tree, parallel=getOption("mc.cores")) # default is permutation=999
signif.full #model is significant## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = samples ~ Land_Coverage_Category + Burn_Severity, data = enviro)
## Df Variance F Pr(>F)
## Model 4 1767.19 7.9608 0.001 ***
## Residual 14 776.95
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#check significance of individual PCAs
signif.axis <- anova.cca(rda_tree, by="axis", parallel=getOption("mc.cores"))
signif.axis #only RDA1 and RDA 2 are sig## Permutation test for rda under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = samples ~ Land_Coverage_Category + Burn_Severity, data = enviro)
## Df Variance F Pr(>F)
## RDA1 1 1300.26 23.4295 0.001 ***
## RDA2 1 399.27 7.1944 0.014 *
## RDA3 1 37.03 0.6672 0.903
## RDA4 1 30.63 0.5519 0.707
## Residual 14 776.95
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#check VIF
vif.cca(rda_tree)## Land_Coverage_CategorySagebrush shrubland
## 1.118041
## Burn_SeverityLow
## 1.846364
## Burn_SeverityModerate
## 1.836781
## Burn_SeverityUnburned
## 1.558868#compare the model to chance with permutation test
anova(rda_tree)## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = samples ~ Land_Coverage_Category + Burn_Severity, data = enviro)
## Df Variance F Pr(>F)
## Model 4 1767.19 7.9608 0.001 ***
## Residual 14 776.95
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#plot again
plot(rda_tree, scaling=2) #this is significant and VIF is good; make this one pretty
#How does feedstock and burn severity influence solid (P concentration, NMR, XANES) and leachate (concentration of different phases and molybdate) chemistry
data_rda<- solid %>%
full_join(leachate_mean) ## Joining with `by = join_by(Parent_ID)`#create dataframe for sample variables
samples<-data_rda[ c(6:17, 27:30)]
#create dataframe for burn condition variables
enviro<-data_rda[ c(2:3)]
#make sure both dataframes have equal rows
all.equal(rownames(samples), rownames(enviro))## [1] TRUE#RDA (this anaysis expects a linear response of each sample to the environmental variables)
rda_tree=rda(samples~ . , data=enviro)
rda_tree## Call: rda(formula = samples ~ Land_Coverage_Category + Burn_Severity,
## data = enviro)
##
## Inertia Proportion Rank
## Total 1.882e+07 1.000e+00
## Constrained 7.940e+06 4.219e-01 4
## Unconstrained 1.088e+07 5.781e-01 12
## Inertia is variance
##
## Eigenvalues for constrained axes:
## RDA1 RDA2 RDA3 RDA4
## 7926038 12958 984 81
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
## 10869073 11608 426 214 92 75 48 17
## PC9 PC10 PC11 PC12
## 8 2 1 0RsquareAdj(rda_tree)## $r.squared
## [1] 0.4218584
##
## $adj.r.squared
## [1] 0.2566751summary(eigenvals(rda_tree, model = "constrained"))## Importance of components:
## RDA1 RDA2 RDA3 RDA4
## Eigenvalue 7.926e+06 1.296e+04 9.840e+02 8.051e+01
## Proportion Explained 9.982e-01 1.632e-03 1.239e-04 1.014e-05
## Cumulative Proportion 9.982e-01 9.999e-01 1.000e+00 1.000e+00#see the factor loadings (correlations of each variable with each pca axis)
scores(rda_tree, choices=1:2, display="species", scaling=2)## RDA1 RDA2
## Ortho 0.31375232 -0.18493684
## Mono -0.19173524 0.13469572
## Di -0.04856659 0.03015374
## Pyro -0.07313601 0.02074476
## XANES_Po -0.27017217 0.03851785
## XANES_Pi_Fe -0.14214762 0.01074035
## XANES_Pi_Na -0.31957053 0.07374803
## XANES_Pi_Ca 0.29969687 -0.14665261
## XANES_Pi_Al -0.03611682 -0.01006229
## XANES_Pi_Mg 0.39900322 -0.05409895
## XANES_Pi_K 0.05464259 0.09400254
## Solid_P_g_kg 0.10505008 0.01294289
## Leachate_P_mg_kgchar_unfilt_Mean 62.56938415 1.27600954
## Leachate_P_mg_kgchar_part_Mean 61.92711377 -1.32452727
## Leachate_P_mg_kgchar_filt0.7_Mean 0.64227038 2.60053681
## MBD_P_mg_kgchar_Mean 0.44612320 1.55961875
## attr(,"const")
## [1] 135.6696#plot it
plot(rda_tree, type='n', scaling=2)
orditorp(rda_tree, display='sp', cex=0.5, scaling=1, col='blue')
text(rda_tree, display='cn', col='red')
#check correlation
pairs.panels(enviro, scale=T)
#can visualize RDA importance
screeplot(rda_tree)
#check significance of all PCAs
signif.full <- anova.cca(rda_tree, parallel=getOption("mc.cores")) # default is permutation=999
signif.full #model is significant## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = samples ~ Land_Coverage_Category + Burn_Severity, data = enviro)
## Df Variance F Pr(>F)
## Model 4 7940061 2.5539 0.059 .
## Residual 14 10881564
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#check significance of individual PCAs
signif.axis <- anova.cca(rda_tree, by="axis", parallel=getOption("mc.cores"))
signif.axis #only RDA1 and RDA 2 are sig## Permutation test for rda under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = samples ~ Land_Coverage_Category + Burn_Severity, data = enviro)
## Df Variance F Pr(>F)
## RDA1 1 7926038 10.1975 0.058 .
## RDA2 1 12958 0.0167 0.999
## RDA3 1 984 0.0013 1.000
## RDA4 1 81 0.0001 0.999
## Residual 14 10881564
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#check VIF
vif.cca(rda_tree)## Land_Coverage_CategorySagebrush shrubland
## 1.118041
## Burn_SeverityLow
## 1.846364
## Burn_SeverityModerate
## 1.836781
## Burn_SeverityUnburned
## 1.558868#compare the model to chance with permutation test
anova(rda_tree)## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = samples ~ Land_Coverage_Category + Burn_Severity, data = enviro)
## Df Variance F Pr(>F)
## Model 4 7940061 2.5539 0.057 .
## Residual 14 10881564
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#plot again
plot(rda_tree, scaling=2) #this is isn't as good with the leachates included
#make a pretty plot
#use
#USE THIS ONE AS AN EXAMPLE!! (has horizons as difference shapes but need to figure out legend)
#levels(enviro$Site) <- c("GR2","GR3","GR5","BAR","SJER","SPR", "PR", "SH")
#eco <- enviro$Site
#bgg <- c("#000080", "#8B0000","#0E5E0B","#606060","#5B94BA","#DF1111", "#069710", "#C9D4D1") # 8 colors for our sites
#shapes <- c(18, 16, 15,17)
#sh <- as.factor(enviro$Horizon)
#plot(rda_tree, type="n", scaling=2)
#points(rda_tree, display="species", pch=20, cex=0.7, col="gray32", scaling=2) #the chemical species from NMR and XANES
#points(rda_tree, display="sites", pch=shapes[sh], cex=1.3, scaling=2, col=bgg[eco]) # the soil samples
#text(rda_tree, scaling=2, display="bp", col="#0868ac", cex=1) # the predictors
#legend("bottomright", legend=levels(eco), bty="n", pch=21, cex=1, col=bg)
#legend("bottomright", legend=c("A", "B", "C", "GR2", "GR3", "GR5", "BAR", "SJER", "SPR", "PR", "SH"), col=c(rep("black",3), "#000080", "#8B0000","#0E5E0B","#606060","#5B94BA","#DF1111", "#069710", "#C9D4D1"), pch=c(16, 15, 17, 16, 16, 16, 16, 16, 16, 16, 16))
#surface horizons only
#Master_CZO_WM_Soil_5.8.2020 <- read.csv("~/Documents/Research Notes/Dissertation General/Master_CZO_WM_Soil_5.8.2020.csv")
#Master_CZO_WM_Soil_5.8.2020$Climo<-factor(Master_CZO_WM_Soil_5.8.2020$Climo, levels=c("WM", "CZO"))
#xas.percent<-subset(Master_CZO_WM_Soil_5.8.2020, Sample.Type=="Soil")
#xas.percent<-xas.percent[ , c("Horizon.Category.C", "Climo", "Site", "EEMT", "AI", "MAT", "MAP","pH.CaCl2", "Clay.Percent", "Fe.g.kgsoil.PMnorm", "Ca.g.kgsoil.PMnorm", "Al.g.kgsoil.PMnorm", "C.g.kgsoil", "P.g.kgsoil.PMnorm", "N.g.kgsoil", "Ca.Pi.Percent", "Al.Pi.Percent", "Fe.Pi.Percent", "XANES.Total.Monoester.Percent", "XANES.Total.Diester.Percent")]
#xas.percent<-subset(xas.percent, Ca.Pi.Percent!="NA")
#colnames(xas.percent)[1:20]<-c("Horizon", "Gradient", "Site","EEMT", "Humidity", "MAT", "MAP","pH", "Clay", "Fe", "Ca", "Al", "C", "P", "N", "Ca-Pi", "Al-Pi", "Fe-Pi", "Monoester", "Diester")
#xas.percent<-subset(xas.percent, Horizon=="A")####stacked spectra and pie charts
#example stacked spectra
NMR_spectra_df_subset <-NMR_spectra %>%
select("Chemical_Shift_parts_per_million", "BSLE_0013-solid", "BSLE_0007-solid","BSLE_0002-solid", "BSLE_0050-solid") %>%
dplyr::rename('solid_13' = 'BSLE_0013-solid',
'solid_7' = 'BSLE_0007-solid',
'solid_2' = 'BSLE_0002-solid',
'solid_50' = 'BSLE_0050-solid') %>%
mutate_all(~replace(., . == 0, NA)) %>%
mutate(solid_13 = solid_13 + 1,
solid_7 = solid_7 + 75,
solid_2 = solid_2 + 150,
solid_50 = solid_50 + 225)
NMR_spectra_df_subset<-melt(NMR_spectra_df_subset, id=c("Chemical_Shift_parts_per_million"))
NMR_spectra_df_subset<-subset(NMR_spectra_df_subset, value!="")
NMR_stacked_spectra_df <- ggplot(NMR_spectra_df_subset, aes(x=Chemical_Shift_parts_per_million, y = value, color = as.factor(variable))) + geom_path()+theme_classic() + ylim(0,300) + scale_x_reverse() +xlim(10,-10) +scale_color_manual(values=c('Black','Black', 'Black', 'Black')) + theme(axis.ticks.y = element_blank(), axis.text.y = element_blank(), legend.position = "none") + ylab("") + xlab("Chemical Shift (ppm)")## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.#cowplot::save_plot("../figures/NMR_stacked_spectra_df.pdf", NMR_stacked_spectra_df, ncol = 1, nrow = 1, base_aspect_ratio= 1:1.7, dpi=72)
pdf(file = "../figures/NMR_stacked_spectra_df.pdf", width = 3.5, height = 7)
print(NMR_stacked_spectra_df)## Warning: Removed 34204 rows containing missing values (`geom_path()`).dev.off()## quartz_off_screen
## 2NMR_spectra_sb_subset <-NMR_spectra %>%
select("Chemical_Shift_parts_per_million", "BSLE_0014-solid", "BSLE_0011-solid","BSLE_0072-solid") %>%
dplyr::rename('solid_14' = 'BSLE_0014-solid',
'solid_11' = 'BSLE_0011-solid',
'solid_72' = 'BSLE_0072-solid') %>%
mutate_all(~replace(., . == 0, NA)) %>%
mutate(solid_14 = solid_14 + 1,
solid_11 = solid_11 + 75,
solid_72 = solid_72 + 150)
NMR_spectra_sb_subset<-melt(NMR_spectra_sb_subset, id=c("Chemical_Shift_parts_per_million"))
NMR_spectra_sb_subset<-subset(NMR_spectra_sb_subset, value!="")
NMR_stacked_spectra_sb <- ggplot(NMR_spectra_sb_subset, aes(x=Chemical_Shift_parts_per_million, y = value, color = as.factor(variable))) + geom_path()+theme_classic() + ylim(0,300) + scale_x_reverse() +xlim(10,-10) +scale_color_manual(values=c('Black','Black', 'Black', 'Black')) + theme(axis.ticks.y = element_blank(), axis.text.y = element_blank(), legend.position = "none") + ylab("") + xlab("Chemical Shift (ppm)")## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.#cowplot::save_plot("../figures/NMR_stacked_spectra_sb.pdf", NMR_stacked_spectra_sb, ncol = 1, nrow = 1, base_aspect_ratio= 1:1.7, dpi=72)
pdf(file = "../figures/NMR_stacked_spectra_sb.pdf", width = 3.5, height = 7)
print(NMR_stacked_spectra_sb)## Warning: Removed 25652 rows containing missing values (`geom_path()`).dev.off()## quartz_off_screen
## 2#pie charts
NMR <- summary_solids_conc %>%
select("Land_Coverage_Category", "Burn_Severity","Ortho_Mean", "Mono_Mean", "Di_Mean", "Pyro_Mean") %>%
dplyr::rename(Ortho = Ortho_Mean,
Mono = Mono_Mean,
Di = Di_Mean,
Pyro = Pyro_Mean)
NMR_melt<-melt(NMR, id=c("Land_Coverage_Category", "Burn_Severity"))
NMR_melt$variable=factor(NMR_melt$variable, levels= c("Ortho", "Pyro", "Mono", "Di"))
DF.raw<-subset(NMR_melt, Land_Coverage_Category=="Douglas-fir forest" & Burn_Severity=="Unburned")
nmr.legend<-ggplot(DF.raw, aes(x = "", y = value, fill = variable)) +
geom_col(color = "black") +
coord_polar(theta = "y") + scale_fill_manual(values = c("#780000", "#c1121f", "#003049", "#669BBC"))+theme_void()+
theme(legend.key.size = unit(1, 'cm'), #change legend key size
legend.key.height = unit(1, 'cm'), legend.title = element_text(size=25),
legend.key.width = unit(1, 'cm'), #change legend key width
legend.text = element_text(size=20))+labs(fill="")
pdf(file = "../figures/nmr.legend.pdf", width = 7, height = 3.5)
print(nmr.legend)
dev.off()## quartz_off_screen
## 2DF.raw.nmr<-ggplot(DF.raw, aes(x = "", y = value, fill = variable)) +
geom_col(color = "black") +
coord_polar(theta = "y") + scale_fill_manual(values = c("#780000", "#c1121f", "#003049", "#669BBC"))+theme_void()+
theme(legend.position = "none")
pdf(file = "../figures/DF.raw.nmr.pdf", width = 1, height = 1)
print(DF.raw.nmr)
dev.off()## quartz_off_screen
## 2DF.low<-subset(NMR_melt, Land_Coverage_Category=="Douglas-fir forest" & Burn_Severity=="Low")
DF.low.nmr<-ggplot(DF.low, aes(x = "", y = value, fill = variable)) +
geom_col(color = "black") +
coord_polar(theta = "y") + scale_fill_manual(values = c("#780000", "#c1121f", "#003049", "#669BBC"))+theme_void()+
theme(legend.position = "none")
pdf(file = "../figures/DF.low.nmr.pdf", width = 1, height = 1)
print(DF.low.nmr)
dev.off()## quartz_off_screen
## 2DF.mod<-subset(NMR_melt, Land_Coverage_Category=="Douglas-fir forest" & Burn_Severity=="Moderate")
DF.mod.nmr<-ggplot(DF.mod, aes(x = "", y = value, fill = variable)) +
geom_col(color = "black") +
coord_polar(theta = "y") + scale_fill_manual(values = c("#780000", "#c1121f", "#003049", "#669BBC"))+theme_void()+
theme(legend.position = "none")
pdf(file = "../figures/DF.mod.nmr.pdf", width = 1, height = 1)
print(DF.mod.nmr)
dev.off()## quartz_off_screen
## 2DF.high<-subset(NMR_melt, Land_Coverage_Category=="Douglas-fir forest" & Burn_Severity=="High")
DF.high.nmr<-ggplot(DF.high, aes(x = "", y = value, fill = variable)) +
geom_col(color = "black") +
coord_polar(theta = "y") + scale_fill_manual(values = c("#780000", "#c1121f", "#003049", "#669BBC"))+theme_void()+
theme(legend.position = "none")
pdf(file = "../figures/DF.high.nmr.pdf", width = 1, height = 1)
print(DF.high.nmr)
dev.off()## quartz_off_screen
## 2SB.raw<-subset(NMR_melt, Land_Coverage_Category=="Sagebrush shrubland" & Burn_Severity=="Unburned")
SB.raw.nmr<-ggplot(SB.raw, aes(x = "", y = value, fill = variable)) +
geom_col(color = "black") +
coord_polar(theta = "y") + scale_fill_manual(values = c("#780000", "#c1121f", "#003049", "#669BBC"))+theme_void()+
theme(legend.position = "none")
pdf(file = "../figures/SB.raw.nmr.pdf", width = 1, height = 1)
print(SB.raw.nmr)
dev.off()## quartz_off_screen
## 2SB.low<-subset(NMR_melt, Land_Coverage_Category=="Sagebrush shrubland" & Burn_Severity=="Low")
SB.low.nmr<-ggplot(SB.low, aes(x = "", y = value, fill = variable)) +
geom_col(color = "black") +
coord_polar(theta = "y") + scale_fill_manual(values = c("#780000", "#c1121f", "#003049", "#669BBC"))+theme_void()+
theme(legend.position = "none")
pdf(file = "../figures/SB.low.nmr.pdf", width = 1, height = 1)
print(SB.low.nmr)
dev.off()## quartz_off_screen
## 2SB.mod<-subset(NMR_melt, Land_Coverage_Category=="Sagebrush shrubland" & Burn_Severity=="Moderate")
SB.mod.nmr<-ggplot(SB.mod, aes(x = "", y = value, fill = variable)) +
geom_col(color = "black") +
coord_polar(theta = "y") + scale_fill_manual(values = c("#780000", "#c1121f", "#003049", "#669BBC"))+theme_void()+
theme(legend.position = "none")
pdf(file = "../figures/SB.mod.nmr.pdf", width = 1, height = 1)
print(SB.mod.nmr)
dev.off()## quartz_off_screen
## 2#example stacked spectra
XANES_spectra_df_subset <-XANES_spectra %>%
select("eV", "BSLE_0013-solid_P-XANES.xmu.nor", "BSLE_0007-solid_P-XANES.xmu.nor","BSLE_0002-solid_P-XANES.xmu.nor", "BSLE_0050-solid_P-XANES.xmu.nor") %>%
dplyr::rename('solid_13' = 'BSLE_0013-solid_P-XANES.xmu.nor',
'solid_7' = 'BSLE_0007-solid_P-XANES.xmu.nor',
'solid_2' = 'BSLE_0002-solid_P-XANES.xmu.nor',
'solid_50' = 'BSLE_0050-solid_P-XANES.xmu.nor') %>%
mutate_all(~replace(., . == 0, NA)) %>%
mutate(solid_13 = solid_13 + 1,
solid_7 = solid_7 + 8,
solid_2 = solid_2 + 16,
solid_50 = solid_50 + 24)
XANES_spectra_df_subset<-melt(XANES_spectra_df_subset, id=c("eV"))
XANES_spectra_df_subset<-subset(XANES_spectra_df_subset, value!="")
XANES_stacked_spectra_df <- ggplot(XANES_spectra_df_subset, aes(x=eV, y = value, color = as.factor(variable))) + geom_path()+theme_classic() + scale_color_manual(values=c('Black','Black', 'Black', 'Black')) + theme(axis.ticks.y = element_blank(), axis.text.y = element_blank(), legend.position = "none") + ylab("") + xlab("Energy (eV)") + xlim(2140, 2200) + ylab ("Normalized Absorption (arb. units)")
pdf(file = "../figures/XANES_stacked_spectra_df.pdf", width = 3, height = 7)
print(XANES_stacked_spectra_df)## Warning: Removed 188 rows containing missing values (`geom_path()`).dev.off()## quartz_off_screen
## 2XANES_spectra_sb_subset <-XANES_spectra %>%
select("eV", "BSLE_0014-solid_P-XANES.xmu.nor", "BSLE_0011-solid_P-XANES.xmu.nor","BSLE_0072-solid_P-XANES.xmu.nor") %>%
dplyr::rename('solid_14' = 'BSLE_0014-solid_P-XANES.xmu.nor',
'solid_11' = 'BSLE_0011-solid_P-XANES.xmu.nor',
'solid_72' = 'BSLE_0072-solid_P-XANES.xmu.nor') %>%
mutate_all(~replace(., . == 0, NA)) %>%
mutate(blank = solid_72) %>%
mutate(solid_14 = solid_14 + 1,
solid_11 = solid_11 + 8,
solid_72 = solid_72 + 16,
blank = blank + 24)
XANES_spectra_sb_subset<-melt(XANES_spectra_sb_subset, id=c("eV"))
XANES_spectra_sb_subset<-subset(XANES_spectra_sb_subset, value!="")
XANES_stacked_spectra_sb <- ggplot(XANES_spectra_sb_subset, aes(x=eV, y = value, color = as.factor(variable))) + geom_path()+theme_classic() + scale_color_manual(values=c('Black','Black', 'Black', 'White')) + theme(axis.ticks.y = element_blank(), axis.text.y = element_blank(), legend.position = "none") + ylab("") + xlab("Energy (eV)") + xlim(2140, 2200) + ylab ("Normalized Absorption (arb. units)")
#cowplot::save_plot("../figures/NMR_stacked_spectra_sb.pdf", NMR_stacked_spectra_sb, ncol = 1, nrow = 1, base_aspect_ratio= 1:1.7, dpi=72)
pdf(file = "../figures/XANES_stacked_spectra_sb.pdf", width = 3, height = 7)
print(XANES_stacked_spectra_sb)## Warning: Removed 188 rows containing missing values (`geom_path()`).dev.off()## quartz_off_screen
## 2#pie charts with Po grouped and inorganic separate metals
XANES <- summary_solids_conc %>%
select("Land_Coverage_Category", "Burn_Severity","XANES_Pi_K_Mean", "XANES_Pi_Mg_Mean", "XANES_Pi_Al_Mean", "XANES_Pi_Ca_Mean", "XANES_Pi_Na_Mean", "XANES_Pi_Fe_Mean", "XANES_Po_Mean") %>%
dplyr::rename(K = XANES_Pi_K_Mean,
Mg = XANES_Pi_Mg_Mean,
Al = XANES_Pi_Al_Mean,
Ca = XANES_Pi_Ca_Mean,
Na = XANES_Pi_Na_Mean,
Fe = XANES_Pi_Fe_Mean,
Po = XANES_Po_Mean)
XANES_melt<-melt(XANES, id=c("Land_Coverage_Category", "Burn_Severity"))
#XANES_melt$variable=factor(XANES_melt$variable, levels= c("Ortho", "Pyro", "Mono", "Di"))
DF.raw<-subset(XANES_melt, Land_Coverage_Category=="Douglas-fir forest" & Burn_Severity=="Unburned")
xanes.legend<-ggplot(DF.raw, aes(x = "", y = value, fill = variable)) +
geom_col(color = "black") +
coord_polar(theta = "y") + scale_fill_manual(values = c("#424a26", "#9B9B7A", "#D9AE94", "#F1DCA7", "#FFCB69", "#D08C60", "#997B66"))+theme_void()+
theme(legend.key.size = unit(1, 'cm'), #change legend key size
legend.key.height = unit(1, 'cm'), legend.title = element_text(size=25),
legend.key.width = unit(1, 'cm'), #change legend key width
legend.text = element_text(size=20))+labs(fill="")
pdf(file = "../figures/xanes.legend.pdf", width = 7, height = 3.5)
print(xanes.legend)
dev.off()## quartz_off_screen
## 2DF.raw.xanes<-ggplot(DF.raw, aes(x = "", y = value, fill = variable)) +
geom_col(color = "black") +
coord_polar(theta = "y") + scale_fill_manual(values = c("#424a26", "#9B9B7A", "#D9AE94", "#F1DCA7", "#FFCB69", "#D08C60", "#997B66"))+theme_void()+
theme(legend.position = "none")
pdf(file = "../figures/DF.raw.xanes.pdf", width = 1, height = 1)
print(DF.raw.xanes)
dev.off()## quartz_off_screen
## 2DF.low<-subset(XANES_melt, Land_Coverage_Category=="Douglas-fir forest" & Burn_Severity=="Low")
DF.low.xanes<-ggplot(DF.low, aes(x = "", y = value, fill = variable)) +
geom_col(color = "black") +
coord_polar(theta = "y") + scale_fill_manual(values = c("#424a26", "#9B9B7A", "#D9AE94", "#F1DCA7", "#FFCB69", "#D08C60", "#997B66"))+theme_void()+
theme(legend.position = "none")
pdf(file = "../figures/DF.low.xanes.pdf", width = 1, height = 1)
print(DF.low.xanes)
dev.off()## quartz_off_screen
## 2DF.mod<-subset(XANES_melt, Land_Coverage_Category=="Douglas-fir forest" & Burn_Severity=="Moderate")
DF.mod.xanes<-ggplot(DF.mod, aes(x = "", y = value, fill = variable)) +
geom_col(color = "black") +
coord_polar(theta = "y") + scale_fill_manual(values = c("#424a26", "#9B9B7A", "#D9AE94", "#F1DCA7", "#FFCB69", "#D08C60", "#997B66"))+theme_void()+
theme(legend.position = "none")
pdf(file = "../figures/DF.mod.xanes.pdf", width = 1, height = 1)
print(DF.mod.xanes)
dev.off()## quartz_off_screen
## 2DF.high<-subset(XANES_melt, Land_Coverage_Category=="Douglas-fir forest" & Burn_Severity=="High")
DF.high.xanes<-ggplot(DF.high, aes(x = "", y = value, fill = variable)) +
geom_col(color = "black") +
coord_polar(theta = "y") + scale_fill_manual(values = c("#424a26", "#9B9B7A", "#D9AE94", "#F1DCA7", "#FFCB69", "#D08C60", "#997B66"))+theme_void()+
theme(legend.position = "none")
pdf(file = "../figures/DF.high.xanes.pdf", width = 1, height = 1)
print(DF.high.xanes)
dev.off()## quartz_off_screen
## 2SB.raw<-subset(XANES_melt, Land_Coverage_Category=="Sagebrush shrubland" & Burn_Severity=="Unburned")
SB.raw.xanes<-ggplot(SB.raw, aes(x = "", y = value, fill = variable)) +
geom_col(color = "black") +
coord_polar(theta = "y") + scale_fill_manual(values = c("#424a26", "#9B9B7A", "#D9AE94", "#F1DCA7", "#FFCB69", "#D08C60", "#997B66"))+theme_void()+
theme(legend.position = "none")
pdf(file = "../figures/SB.raw.xanes.pdf", width = 1, height = 1)
print(SB.raw.xanes)
dev.off()## quartz_off_screen
## 2SB.low<-subset(XANES_melt, Land_Coverage_Category=="Sagebrush shrubland" & Burn_Severity=="Low")
SB.low.xanes<-ggplot(SB.low, aes(x = "", y = value, fill = variable)) +
geom_col(color = "black") +
coord_polar(theta = "y") + scale_fill_manual(values = c("#424a26", "#9B9B7A", "#D9AE94", "#F1DCA7", "#FFCB69", "#D08C60", "#997B66"))+theme_void()+
theme(legend.position = "none")
pdf(file = "../figures/SB.low.xanes.pdf", width = 1, height = 1)
print(SB.low.xanes)
dev.off()## quartz_off_screen
## 2SB.mod<-subset(XANES_melt, Land_Coverage_Category=="Sagebrush shrubland" & Burn_Severity=="Moderate")
SB.mod.xanes<-ggplot(SB.mod, aes(x = "", y = value, fill = variable)) +
geom_col(color = "black") +
coord_polar(theta = "y") + scale_fill_manual(values = c("#424a26", "#9B9B7A", "#D9AE94", "#F1DCA7", "#FFCB69", "#D08C60", "#997B66"))+theme_void()+
theme(legend.position = "none")
pdf(file = "../figures/SB.mod.xanes.pdf", width = 1, height = 1)
print(SB.mod.xanes)
dev.off()## quartz_off_screen
## 2#SEM
#reference: https://stats.oarc.ucla.edu/r/seminars/rsem/ #followed Model 4A
#make dataframe with solid samples (leachate avg) to prevent using same number so many times
leachate_mean <- leachate %>%
dplyr::group_by(Parent_ID) %>%
dplyr::summarize(
Leachate_P_mg_kgchar_unfilt_Mean = mean(Leachate_P_mg_kgchar_unfilt, na.rm = TRUE),
Leachate_P_mg_kgchar_part_Mean = mean(Leachate_P_mg_kgchar_part, na.rm = TRUE),
Leachate_P_mg_kgchar_filt0.7_Mean = mean(Leachate_P_mg_kgchar_filt0.7, na.rm = TRUE),
MBD_P_mg_kgchar_Mean = mean(MBD_P_mg_kgchar, na.rm = TRUE),
pH_Mean = mean(pH)
)
solid <- data_solids %>%
select(1,3, 5, 10, 11, 31:34, 88, 96:101, 111:120)
data_rda<- solid %>%
full_join(leachate_mean) %>%
filter(Solid_Na_g_kg != "") #remove NAs## Joining with `by = join_by(Parent_ID)`model <- '
# Measurement Model
f1 =~ Land_Coverage_Category + Char_Max_Temp + Burn_Duration
f2 =~ Burn_Severity
# Structural Model
f1 ~ a*f2
y1 ~ b*f1
# Residual Variances
Land_Coverage_Category ~~ Land_Coverage_Category
Char_Max_Temp ~~ Char_Max_Temp
Burn_Duration ~~ Burn_Duration
Burn_Severity ~~ Burn_Severity
'
Leachate_P_mg_kgchar_unfilt_Mean ~ 1 + Solid_P_g_kg + Mono + Di + XANES_Pi_Ca + XANES_Pi_Mg + Burn_Severity + Land_Coverage_Category + Burn_Duration + Char_Max_Temp## Leachate_P_mg_kgchar_unfilt_Mean ~ 1 + Solid_P_g_kg + Mono +
## Di + XANES_Pi_Ca + XANES_Pi_Mg + Burn_Severity + Land_Coverage_Category +
## Burn_Duration + Char_Max_TempSolid_P_g_kg ~ 1 + Burn_Severity + Land_Coverage_Category + Burn_Duration + Char_Max_Temp## Solid_P_g_kg ~ 1 + Burn_Severity + Land_Coverage_Category + Burn_Duration +
## Char_Max_TempMono ~ 1 + Burn_Severity + Land_Coverage_Category + Burn_Duration + Char_Max_Temp## Mono ~ 1 + Burn_Severity + Land_Coverage_Category + Burn_Duration +
## Char_Max_TempDi ~ 1 + Burn_Severity + Land_Coverage_Category + Burn_Duration + Char_Max_Temp## Di ~ 1 + Burn_Severity + Land_Coverage_Category + Burn_Duration +
## Char_Max_TempXANES_Pi_Ca ~ 1 + Burn_Severity + Land_Coverage_Category + Burn_Duration + Char_Max_Temp## XANES_Pi_Ca ~ 1 + Burn_Severity + Land_Coverage_Category + Burn_Duration +
## Char_Max_TempXANES_Pi_Mg ~ 1 + Burn_Severity + Land_Coverage_Category + Burn_Duration + Char_Max_Temp## XANES_Pi_Mg ~ 1 + Burn_Severity + Land_Coverage_Category + Burn_Duration +
## Char_Max_Tempmodel<-'
Solid_P_g_kg ~ 1 + Burn_Severity
Mono ~ 1 + Burn_Severity
Di ~ 1 + Burn_Severity
XANES_Pi_Ca ~ 1 + Burn_Severity
XANES_Pi_Mg ~ 1 + Burn_Severity
Burn_Severity ~ 1 + Land_Coverage_Category + Burn_Duration + Char_Max_Temp
Char_Max_Temp ~ Land_Coverage_Category + Burn_Duration
'
fit <- sem(model, data = data_rda)## Warning in lav_data_full(data = data, group = group, cluster = cluster, :
## lavaan WARNING: some observed variances are (at least) a factor 1000 times
## larger than others; use varTable(fit) to investigatesummary(fit, fit.measures=TRUE)## lavaan 0.6.16 ended normally after 181 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 34
##
## Number of observations 18
##
## Model Test User Model:
##
## Test statistic 63.543
## Degrees of freedom 15
## P-value (Chi-square) 0.000
##
## Model Test Baseline Model:
##
## Test statistic 185.488
## Degrees of freedom 35
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.677
## Tucker-Lewis Index (TLI) 0.247
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -418.634
## Loglikelihood unrestricted model (H1) NA
##
## Akaike (AIC) 905.268
## Bayesian (BIC) 935.541
## Sample-size adjusted Bayesian (SABIC) 831.069
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.424
## 90 Percent confidence interval - lower 0.320
## 90 Percent confidence interval - upper 0.534
## P-value H_0: RMSEA <= 0.050 0.000
## P-value H_0: RMSEA >= 0.080 1.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.158
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Regressions:
## Estimate Std.Err z-value P(>|z|)
## Solid_P_g_kg ~
## Burn_Severity -0.684 0.968 -0.707 0.480
## Mono ~
## Burn_Severity 7.760 2.295 3.381 0.001
## Di ~
## Burn_Severity 2.949 0.715 4.126 0.000
## XANES_Pi_Ca ~
## Burn_Severity -11.707 2.779 -4.212 0.000
## XANES_Pi_Mg ~
## Burn_Severity -4.667 4.299 -1.086 0.278
## Burn_Severity ~
## Lnd_Cvrg_Ctgry 0.042 0.412 0.103 0.918
## Burn_Duration -0.001 0.001 -1.077 0.282
## Char_Max_Temp -0.002 0.001 -2.095 0.036
## Char_Max_Temp ~
## Lnd_Cvrg_Ctgry 29.882 99.784 0.299 0.765
## Burn_Duration 0.501 0.127 3.941 0.000
##
## Covariances:
## Estimate Std.Err z-value P(>|z|)
## .Solid_P_g_kg ~~
## .Mono -21.931 10.877 -2.016 0.044
## .Di -4.552 3.167 -1.437 0.151
## .XANES_Pi_Ca 21.671 12.665 1.711 0.087
## .XANES_Pi_Mg 54.259 22.020 2.464 0.014
## .Mono ~~
## .Di 23.678 9.005 2.629 0.009
## .XANES_Pi_Ca -92.587 35.092 -2.638 0.008
## .XANES_Pi_Mg -95.291 48.076 -1.982 0.047
## .Di ~~
## .XANES_Pi_Ca -27.776 10.774 -2.578 0.010
## .XANES_Pi_Mg -25.395 14.526 -1.748 0.080
## .XANES_Pi_Ca ~~
## .XANES_Pi_Mg 122.063 58.965 2.070 0.038
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|)
## .Solid_P_g_kg 6.107 2.509 2.434 0.015
## .Mono -8.444 5.951 -1.419 0.156
## .Di -4.745 1.853 -2.561 0.010
## .XANES_Pi_Ca 64.778 7.206 8.990 0.000
## .XANES_Pi_Mg 31.576 11.146 2.833 0.005
## .Burn_Severity 3.592 0.698 5.146 0.000
## .Char_Max_Temp 254.362 158.596 1.604 0.109
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .Solid_P_g_kg 17.123 5.708 3.000 0.003
## .Mono 96.287 32.096 3.000 0.003
## .Di 9.336 3.112 3.000 0.003
## .XANES_Pi_Ca 141.180 47.060 3.000 0.003
## .XANES_Pi_Mg 337.761 112.587 3.000 0.003
## .Burn_Severity 0.502 0.167 3.000 0.003
## .Char_Max_Temp 29645.894 9881.965 3.000 0.003modindices(fit, sort=TRUE)## lhs op rhs mi epc sepc.lv
## 94 Char_Max_Temp ~ XANES_Pi_Ca 10.746 11.939 11.939
## 93 Char_Max_Temp ~ Di 9.911 -44.588 -44.588
## 92 Char_Max_Temp ~ Mono 9.856 -13.846 -13.846
## 87 Burn_Severity ~ Mono 8.165 -0.068 -0.068
## 97 Land_Coverage_Category ~ Solid_P_g_kg 7.567 0.064 0.064
## 109 Burn_Duration ~ XANES_Pi_Mg 7.452 13.558 13.558
## 89 Burn_Severity ~ XANES_Pi_Ca 7.164 0.053 0.053
## 101 Land_Coverage_Category ~ XANES_Pi_Mg 5.624 0.012 0.012
## 85 XANES_Pi_Mg ~ Burn_Duration 5.609 0.024 0.024
## 90 Burn_Severity ~ XANES_Pi_Mg 5.516 0.030 0.030
## 56 Solid_P_g_kg ~ Land_Coverage_Category 5.464 3.443 3.443
## 88 Burn_Severity ~ Di 5.337 -0.178 -0.178
## 40 Solid_P_g_kg ~~ Burn_Severity 3.020 -1.106 -1.106
## 48 XANES_Pi_Mg ~~ Burn_Severity 2.986 4.796 4.796
## 57 Solid_P_g_kg ~ Burn_Duration 2.968 -0.004 -0.004
## 105 Burn_Duration ~ Solid_P_g_kg 2.585 35.466 35.466
## 63 Mono ~ Land_Coverage_Category 2.413 4.119 4.119
## 42 Mono ~~ Burn_Severity 2.370 -1.764 -1.764
## 62 Mono ~ Char_Max_Temp 1.955 -0.009 -0.009
## 55 Solid_P_g_kg ~ Char_Max_Temp 1.681 -0.005 -0.005
## 100 Land_Coverage_Category ~ XANES_Pi_Ca 1.595 0.010 0.010
## 64 Mono ~ Burn_Duration 1.315 -0.005 -0.005
## 108 Burn_Duration ~ XANES_Pi_Ca 1.223 8.496 8.496
## 106 Burn_Duration ~ Mono 1.149 -9.970 -9.970
## 83 XANES_Pi_Mg ~ Char_Max_Temp 1.010 0.016 0.016
## 77 XANES_Pi_Ca ~ Land_Coverage_Category 0.897 3.312 3.312
## 95 Char_Max_Temp ~ XANES_Pi_Mg 0.840 2.158 2.158
## 91 Char_Max_Temp ~ Solid_P_g_kg 0.824 9.496 9.496
## 86 Burn_Severity ~ Solid_P_g_kg 0.754 0.049 0.049
## 71 Di ~ Burn_Duration 0.673 0.001 0.001
## 76 XANES_Pi_Ca ~ Char_Max_Temp 0.537 0.006 0.006
## 98 Land_Coverage_Category ~ Mono 0.500 -0.007 -0.007
## 99 Land_Coverage_Category ~ Di 0.404 -0.020 -0.020
## 78 XANES_Pi_Ca ~ Burn_Duration 0.341 -0.003 -0.003
## 45 Di ~~ Char_Max_Temp 0.340 -43.373 -43.373
## 107 Burn_Duration ~ Di 0.297 -16.287 -16.287
## 44 Di ~~ Burn_Severity 0.174 0.169 0.169
## 84 XANES_Pi_Mg ~ Land_Coverage_Category 0.109 -2.117 -2.117
## 70 Di ~ Land_Coverage_Category 0.109 -0.308 -0.308
## 46 XANES_Pi_Ca ~~ Burn_Severity 0.059 0.368 0.368
## 69 Di ~ Char_Max_Temp 0.007 0.000 0.000
## sepc.all sepc.nox
## 94 0.824 0.824
## 93 -0.783 -0.783
## 92 -0.716 -0.716
## 87 -0.852 -0.852
## 97 0.597 0.597
## 109 0.732 0.732
## 89 0.880 0.880
## 101 0.524 0.524
## 85 0.443 0.001
## 90 0.566 0.566
## 56 0.368 0.821
## 88 -0.752 -0.752
## 40 -0.377 -0.377
## 48 0.368 0.368
## 57 -0.334 -0.001
## 105 0.423 0.423
## 63 0.147 0.328
## 42 -0.254 -0.254
## 62 -0.175 -0.175
## 55 -0.270 -0.270
## 100 0.381 0.381
## 64 -0.134 0.000
## 108 0.405 0.405
## 106 -0.356 -0.356
## 83 0.202 0.202
## 77 0.089 0.198
## 95 0.169 0.169
## 91 0.164 0.164
## 86 0.205 0.205
## 71 0.099 0.000
## 76 0.091 0.091
## 98 -0.194 -0.194
## 99 -0.190 -0.190
## 78 -0.067 0.000
## 45 -0.082 -0.082
## 107 -0.197 -0.197
## 44 0.078 0.078
## 84 -0.050 -0.112
## 70 -0.032 -0.072
## 46 0.044 0.044
## 69 0.011 0.011model<-'
Leachate_P_mg_kgchar_unfilt_Mean ~ 1 + Solid_P_g_kg + Mono + Di + XANES_Pi_Ca + XANES_Pi_Mg + Burn_Severity
Solid_P_g_kg ~ 1 + Burn_Severity
Mono ~ 1 + Burn_Severity
Di ~ 1 + Burn_Severity
XANES_Pi_Ca ~ 1 + Burn_Severity
XANES_Pi_Mg ~ 1 + Burn_Severity
'
fit <- sem(model, data = data_rda)## Warning in lav_data_full(data = data, group = group, cluster = cluster, :
## lavaan WARNING: some observed variances are (at least) a factor 1000 times
## larger than others; use varTable(fit) to investigate## Warning in lav_data_full(data = data, group = group, cluster = cluster, : lavaan WARNING: some observed variances are larger than 1000000
## lavaan NOTE: use varTable(fit) to investigatesummary(fit, fit.measures=TRUE)## lavaan 0.6.16 ended normally after 1 iteration
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 23
##
## Number of observations 18
##
## Model Test User Model:
##
## Test statistic 60.715
## Degrees of freedom 10
## P-value (Chi-square) 0.000
##
## Model Test Baseline Model:
##
## Test statistic 154.464
## Degrees of freedom 21
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.620
## Tucker-Lewis Index (TLI) 0.202
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -464.191
## Loglikelihood unrestricted model (H1) NA
##
## Akaike (AIC) 974.382
## Bayesian (BIC) 994.861
## Sample-size adjusted Bayesian (SABIC) 924.189
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.531
## 90 Percent confidence interval - lower 0.407
## 90 Percent confidence interval - upper 0.663
## P-value H_0: RMSEA <= 0.050 0.000
## P-value H_0: RMSEA >= 0.080 1.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.267
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Regressions:
## Estimate Std.Err z-value P(>|z|)
## Leachate_P_mg_kgchar_unfilt_Mean ~
## Solid_P_g_kg 748.443 66.843 11.197 0.000
## Mono 18.575 28.188 0.659 0.510
## Di -98.463 90.523 -1.088 0.277
## XANES_Pi_Ca -30.447 23.279 -1.308 0.191
## XANES_Pi_Mg -13.373 15.050 -0.889 0.374
## Burn_Severity 11.290 525.112 0.021 0.983
## Solid_P_g_kg ~
## Burn_Severity -0.684 0.968 -0.707 0.480
## Mono ~
## Burn_Severity 7.760 2.295 3.381 0.001
## Di ~
## Burn_Severity 2.949 0.715 4.126 0.000
## XANES_Pi_Ca ~
## Burn_Severity -11.707 2.779 -4.212 0.000
## XANES_Pi_Mg ~
## Burn_Severity -4.667 4.299 -1.086 0.278
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|)
## .Lcht_P_mg_k__M -213.388 1847.683 -0.115 0.908
## .Solid_P_g_kg 6.107 2.509 2.434 0.015
## .Mono -8.444 5.951 -1.419 0.156
## .Di -4.745 1.853 -2.561 0.010
## .XANES_Pi_Ca 64.778 7.206 8.990 0.000
## .XANES_Pi_Mg 31.576 11.146 2.833 0.005
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .Lcht_P_mg_k__M 1377085.799 459028.600 3.000 0.003
## .Solid_P_g_kg 17.123 5.708 3.000 0.003
## .Mono 96.287 32.096 3.000 0.003
## .Di 9.336 3.112 3.000 0.003
## .XANES_Pi_Ca 141.180 47.060 3.000 0.003
## .XANES_Pi_Mg 337.761 112.587 3.000 0.003modindices(fit, sort=TRUE)## lhs op rhs mi epc sepc.lv
## 36 Mono ~~ XANES_Pi_Ca 11.351 -92.587 -92.587
## 58 XANES_Pi_Ca ~ Mono 11.351 -0.962 -0.962
## 49 Mono ~ XANES_Pi_Ca 11.351 -0.656 -0.656
## 35 Mono ~~ Di 11.226 23.678 23.678
## 53 Di ~ Mono 11.226 0.246 0.246
## 48 Mono ~ Di 11.226 2.536 2.536
## 38 Di ~~ XANES_Pi_Ca 10.536 -27.776 -27.776
## 59 XANES_Pi_Ca ~ Di 10.536 -2.975 -2.975
## 54 Di ~ XANES_Pi_Ca 10.536 -0.197 -0.197
## 45 Solid_P_g_kg ~ XANES_Pi_Mg 9.163 0.161 0.161
## 34 Solid_P_g_kg ~~ XANES_Pi_Mg 9.163 54.259 54.259
## 62 XANES_Pi_Mg ~ Solid_P_g_kg 9.163 3.169 3.169
## 61 XANES_Pi_Mg ~ Leachate_P_mg_kgchar_unfilt_Mean 7.661 0.004 0.004
## 41 Solid_P_g_kg ~ Leachate_P_mg_kgchar_unfilt_Mean 6.033 -0.004 -0.004
## 65 XANES_Pi_Mg ~ XANES_Pi_Ca 5.624 0.865 0.865
## 60 XANES_Pi_Ca ~ XANES_Pi_Mg 5.624 0.361 0.361
## 40 XANES_Pi_Ca ~~ XANES_Pi_Mg 5.624 122.063 122.063
## 47 Mono ~ Solid_P_g_kg 5.251 -1.281 -1.281
## 31 Solid_P_g_kg ~~ Mono 5.251 -21.931 -21.931
## 42 Solid_P_g_kg ~ Mono 5.251 -0.228 -0.228
## 63 XANES_Pi_Mg ~ Mono 5.026 -0.990 -0.990
## 37 Mono ~~ XANES_Pi_Mg 5.026 -95.291 -95.291
## 50 Mono ~ XANES_Pi_Mg 5.026 -0.282 -0.282
## 46 Mono ~ Leachate_P_mg_kgchar_unfilt_Mean 4.065 -0.001 -0.001
## 64 XANES_Pi_Mg ~ Di 3.681 -2.720 -2.720
## 39 Di ~~ XANES_Pi_Mg 3.681 -25.395 -25.395
## 55 Di ~ XANES_Pi_Mg 3.681 -0.075 -0.075
## 33 Solid_P_g_kg ~~ XANES_Pi_Ca 3.497 21.671 21.671
## 44 Solid_P_g_kg ~ XANES_Pi_Ca 3.497 0.153 0.153
## 57 XANES_Pi_Ca ~ Solid_P_g_kg 3.497 1.266 1.266
## 56 XANES_Pi_Ca ~ Leachate_P_mg_kgchar_unfilt_Mean 3.175 0.002 0.002
## 32 Solid_P_g_kg ~~ Di 2.333 -4.552 -4.552
## 43 Solid_P_g_kg ~ Di 2.333 -0.488 -0.488
## 52 Di ~ Solid_P_g_kg 2.333 -0.266 -0.266
## 51 Di ~ Leachate_P_mg_kgchar_unfilt_Mean 0.623 0.000 0.000
## 66 Burn_Severity ~ Leachate_P_mg_kgchar_unfilt_Mean 0.000 0.000 0.000
## 67 Burn_Severity ~ Solid_P_g_kg 0.000 0.000 0.000
## 71 Burn_Severity ~ XANES_Pi_Mg 0.000 0.000 0.000
## 24 Burn_Severity ~~ Burn_Severity 0.000 0.000 0.000
## 70 Burn_Severity ~ XANES_Pi_Ca 0.000 0.000 0.000
## 68 Burn_Severity ~ Mono 0.000 0.000 0.000
## 69 Burn_Severity ~ Di 0.000 0.000 0.000
## sepc.all sepc.nox
## 36 -0.794 -0.794
## 58 -0.721 -0.721
## 49 -0.875 -0.875
## 35 0.790 0.790
## 53 0.724 0.724
## 48 0.862 0.862
## 38 -0.765 -0.765
## 59 -0.757 -0.757
## 54 -0.773 -0.773
## 45 0.726 0.726
## 34 0.713 0.713
## 62 0.701 0.701
## 61 0.678 0.678
## 41 -3.426 -3.426
## 65 0.763 0.763
## 60 0.409 0.409
## 40 0.559 0.559
## 47 -0.428 -0.428
## 31 -0.540 -0.540
## 42 -0.681 -0.681
## 63 -0.655 -0.655
## 37 -0.528 -0.528
## 50 -0.427 -0.427
## 46 -0.398 -0.398
## 64 -0.611 -0.611
## 39 -0.452 -0.452
## 55 -0.335 -0.335
## 33 0.441 0.441
## 44 0.613 0.613
## 57 0.317 0.317
## 56 0.321 0.321
## 32 -0.360 -0.360
## 43 -0.495 -0.495
## 52 -0.262 -0.262
## 51 -0.143 -0.143
## 66 0.000 0.000
## 67 0.000 0.000
## 71 0.000 0.000
## 24 0.000 0.000
## 70 0.000 0.000
## 68 0.000 0.000
## 69 0.000 0.000#reference: https://kevintshoemaker.github.io/NRES-746/SEM.RMarkdown.html
#making a structural regression
model<-'
chem =~ Solid_P_g_kg + Mono + Di + XANES_Pi_Ca + XANES_Pi_Mg
mobilization =~ Leachate_P_mg_kgchar_filt0.7_Mean + MBD_P_mg_kgchar_Mean
mobilization ~ chem + Burn_Severity
'
fit<-sem(model, data = data_rda)## Warning in lav_object_post_check(object): lavaan WARNING: some estimated ov
## variances are negativesummary(fit, fit.measures=TRUE, standardized=TRUE)## lavaan 0.6.16 ended normally after 335 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 16
##
## Number of observations 18
##
## Model Test User Model:
##
## Test statistic 41.992
## Degrees of freedom 19
## P-value (Chi-square) 0.002
##
## Model Test Baseline Model:
##
## Test statistic 156.638
## Degrees of freedom 28
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.821
## Tucker-Lewis Index (TLI) 0.737
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -488.780
## Loglikelihood unrestricted model (H1) NA
##
## Akaike (AIC) 1009.559
## Bayesian (BIC) 1023.805
## Sample-size adjusted Bayesian (SABIC) 974.642
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.259
## 90 Percent confidence interval - lower 0.153
## 90 Percent confidence interval - upper 0.366
## P-value H_0: RMSEA <= 0.050 0.003
## P-value H_0: RMSEA >= 0.080 0.994
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.226
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## chem =~
## Solid_P_g_kg 1.000 2.066 0.493
## Mono -5.738 2.494 -2.301 0.021 -11.856 -0.945
## Di -1.911 0.836 -2.286 0.022 -3.948 -0.926
## XANES_Pi_Ca 7.601 3.311 2.295 0.022 15.706 0.938
## XANES_Pi_Mg 5.381 2.910 1.849 0.064 11.120 0.586
## mobilization =~
## Lcht_P___0.7_M 1.000 98.360 0.829
## MBD_P_mg_kgc_M 1.014 2.044 0.496 0.620 99.733 1.158
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## mobilization ~
## chem 0.662 8.050 0.082 0.934 0.014 0.014
## Burn_Severity -4.959 27.075 -0.183 0.855 -0.050 -0.051
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .Solid_P_g_kg 13.328 4.507 2.957 0.003 13.328 0.757
## .Mono 16.870 9.445 1.786 0.074 16.870 0.107
## .Di 2.581 1.212 2.129 0.033 2.581 0.142
## .XANES_Pi_Ca 33.668 17.467 1.928 0.054 33.668 0.120
## .XANES_Pi_Mg 236.213 80.621 2.930 0.003 236.213 0.656
## .Lcht_P___0.7_M 4417.761 19515.513 0.226 0.821 4417.761 0.313
## .MBD_P_mg_kgc_M -2529.317 20024.632 -0.126 0.899 -2529.317 -0.341
## chem 4.270 3.905 1.093 0.274 1.000 1.000
## .mobilization 9647.787 19723.108 0.489 0.625 0.997 0.997# Define the SEM model
model <- '
# Measurement Model
Chem =~ Solid_P_g_kg + Mono + Di + XANES_Pi_Ca + XANES_Pi_Mg
Mobilization =~ Leachate_P_mg_kgchar_unfilt_Mean + MBD_P_mg_kgchar_Mean
# Structural Model
Mobilization ~ c*Burn_Severity + d*Chem + e*Burn_Severity*Chem
# Residual Variances
Land_Coverage_Category ~~ Land_Coverage_Category
Burn_Duration ~~ Burn_Duration
Char_Max_Temp ~~ Char_Max_Temp
Solid_P_g_kg ~~ Solid_P_g_kg
Mono ~~ Mono
Di ~~ Di
XANES_Pi_Ca ~~ XANES_Pi_Ca
XANES_Pi_Mg ~~ XANES_Pi_Mg
Leachate_P_mg_kgchar_unfilt_Mean ~~ Leachate_P_mg_kgchar_unfilt_Mean
MBD_P_mg_kgchar_Mean ~~ MBD_P_mg_kgchar_Mean
'
# Fit the SEM model
fit <- sem(model, data = data_rda)## Warning in lav_data_full(data = data, group = group, cluster = cluster, : lavaan WARNING: some observed variances are (at least) a factor 1000 times larger than others; use varTable(fit) to investigate
## Warning in lav_data_full(data = data, group = group, cluster = cluster, : lavaan WARNING: some observed variances are larger than 1000000
## lavaan NOTE: use varTable(fit) to investigate## Warning in lav_model_vcov(lavmodel = lavmodel, lavsamplestats = lavsamplestats, : lavaan WARNING:
## Could not compute standard errors! The information matrix could
## not be inverted. This may be a symptom that the model is not
## identified.# Summarize the results
summary(fit, fit.measures=TRUE, standardized=TRUE)## lavaan 0.6.16 ended normally after 119 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 19
##
## Number of observations 18
##
## Model Test User Model:
##
## Test statistic 175.445
## Degrees of freedom 46
## P-value (Chi-square) 0.000
##
## Model Test Baseline Model:
##
## Test statistic 273.519
## Degrees of freedom 55
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.408
## Tucker-Lewis Index (TLI) 0.292
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -834.409
## Loglikelihood unrestricted model (H1) NA
##
## Akaike (AIC) 1706.819
## Bayesian (BIC) 1723.736
## Sample-size adjusted Bayesian (SABIC) 1665.355
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.395
## 90 Percent confidence interval - lower 0.334
## 90 Percent confidence interval - upper 0.458
## P-value H_0: RMSEA <= 0.050 0.000
## P-value H_0: RMSEA >= 0.080 1.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.367
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## Chem =~
## Solid_P_g_kg 1.000 2.188 0.522
## Mono -5.444 NA -11.914 -0.950
## Di -1.793 NA -3.924 -0.921
## XANES_Pi_Ca 7.126 NA 15.596 0.931
## XANES_Pi_Mg 5.188 NA 11.353 0.598
## Mobilization =~
## Lcht_P_mg_k__M 1.000 2123.611 0.619
## MBD_P_mg_kgc_M -0.005 NA -10.906 -0.126
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## Mobilization ~
## Burn_Svrty (c) 1005.436 NA 0.473 0.477
## Chem (d) 852.823 NA 0.879 0.879
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## Lnd_Cvrg_Ctgry 0.201 NA 0.201 1.000
## Burn_Duration 123577.396 NA 123577.396 1.000
## Char_Max_Temp 58868.714 NA 58868.714 1.000
## .Solid_P_g_kg 12.809 NA 12.809 0.728
## .Mono 15.484 NA 15.484 0.098
## .Di 2.767 NA 2.767 0.152
## .XANES_Pi_Ca 37.119 NA 37.119 0.132
## .XANES_Pi_Mg 230.993 NA 230.993 0.642
## .Lcht_P_mg_k__M 7276330.946 NA 7276330.946 0.617
## .MBD_P_mg_kgc_M 7346.177 NA 7346.177 0.984
## Chem 4.789 NA 1.000 1.000
## .Mobilization 0.678 NA 0.000 0.000#haven’t updated below here yet (NMR and XANES code) fix
#spectra stacked #spectra <- read.csv(“~/Library/CloudStorage/OneDrive-SharedLibraries-PNNL/RC-3, River Corridor SFA - Presentations/ESA_2023/Morgan/spectra.csv”)
#filt<-spectra[ -c(6:9)] #filt\(P_BSLE_0007_filt07_sample<-filt\)P_BSLE_0007_filt07_sample+7 #filt\(P_BSLE_0002_filt07_sample<-filt\)P_BSLE_0002_filt07_sample+14 #filt\(P_BSLE_0050_filt07_sample<-filt\)P_BSLE_0050_filt07_sample+21 #filt<-melt(filt, id=c(“energy”))
#ggplot(filt, aes(x=energy, y=value))+geom_line(aes(color = variable), show.legend = FALSE)+theme_classic()+ylab(“Noramlized Absorption (arb. units)”)+ xlab(“Energy (eV)”)+theme(text = element_text(size = 15))+xlim(2140,2200)+scale_color_manual(values=c(‘Black’,‘Black’, ‘Black’, ‘Black’))
#plot spectra of samples #BSLE_solid_sample_spectra_normalizedforfinalLCF_12.1.22 <- read.csv(“~/PNNL/RC-3, River Corridor SFA - Burn Severity Experiments/Burn Severity Experiment/Data/P_all data streams/P XAS/Processed_final/BSLE_Solids_Final_Fits/BSLE_solid_sample_spectra_normalizedforfinalLCF_12.1.22.csv”)
#DFL<-BSLE_solid_sample_spectra_normalizedforfinalLCF_12.1.22[ c(25, 20, 19, 13, 1)] #DFL<-melt(DFL, id=c(“energy”))
#ggplot(DFL, aes(x=energy, y=value))+geom_line()+theme_classic()+ylab(“Noramlized Absorption (arb. units)”)+ xlab(“Energy (eV)”)+facet_grid(variable~.)+theme(text = element_text(size = 15))+xlim(2140,2200)
#try to stack them differently #DFL<-BSLE_solid_sample_spectra_normalizedforfinalLCF_12.1.22[ c(25, 20, 19, 13, 1)] #DFL\(P_BSLE_0073_solidchar_036_fl_4s_avg<-DFL\)P_BSLE_0073_solidchar_036_fl_4s_avg+9 #DFL\(P_BSLE_0052_solidchar_035_fl_6sc_avg<-DFL\)P_BSLE_0052_solidchar_035_fl_6sc_avg+6 #DFL\(P_BSLE_0051_solidchar_034_fl_6s_avg<-DFL\)P_BSLE_0051_solidchar_034_fl_6s_avg+3 #DFL<-melt(DFL, id=c(“energy”))
#DFL.sample.xas.spectra.gradient<-ggplot(DFL, aes(x=energy, y=value))+geom_line(aes(color = variable), show.legend = FALSE)+theme_classic()+ylab(“Noramlized Absorption (arb. units)”)+ xlab(“Energy (eV)”)+theme(text = element_text(size = 15))+xlim(2140,2200)+scale_color_manual(values=c(‘Black’,‘Black’, ‘Black’, ‘Black’))